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Abstract 

Fault detection and isolation (FDI) is a task to deduce 
from observed variable of the system if any component is 
faulty, to locate the faulty components and also to 
estimate the fault magnitude present in the system. This 
paper provides a systematic method of fault diagnosis to 
detect leak in the three-tank process. The proposed 
scheme makes use of structured residual approach for 
detection, isolation and estimation of faults acting on the 
process [1,2,3]. This technique includes residual 
generation and residual evaluation. A literature review 
showed that the conventional fault diagnosis methods like 
the ordinary Chi-square (¥2) test method, generalized 
likelihood ratio test have limitations such as the “false 
alarm” problem. From the results it is inferred that the 
proposed FDI scheme diagnoses better when compared to 
other conventional methods. 

Keywords: Fault detection and Diagnosis, Residual 

Generation, Structured Residual Approach 

1. Introduction 

A fault is any kind of malfunction in a system that is a 
deviation from the normal behavior in the plant or its 
instrumentation [4]. Fault detection is the indication that 
something is going wrong in the monitored system. Fault 
isolation is the determination of the exact location of the 
fault (the component which is faulty) [5,6]. Detection and 
isolation of faults are very important tasks in intelligent 
control systems because a fault can lead to a reduction in 
performance or even to breakdowns or catastrophes. 
Many different approaches to Fault Detection and 
Isolation (FDI) have been proposed over the years. They 
are based either on the model of the system (analytical 
methods) [6, 7, 8, 9] or on the knowledge about the 
system (heuristic methods) [10,11,12, 13, 14]. One way to 
detect faults in a system is by means of analytical 
redundancy [5, 9]. This consists of comparing the 
behaviour of the real system and the behaviour of a model 
of the system. In an ideal case, the system and the model 
behave exactly the same and a fault is detected when the 
behaviours are different, but usually there are differences 
between the behaviours of the system and the model. 



Among the various FDI schemes, the structured 
residual approach (SRA) proposed by Gertler, J., 
Staroswiecki, M. and Shen, M [3] is powerful in 
isolating faults. In this paper the structured residual 
approach is applied to obtain estimates of the 
envelopes for detecting faults in a Three-Tank System 
[6] system which is formed by a sequence of three 
interconnected tanks and has been declared as a 
benchmark. This system has the typical characteristics 
of tanks, pipelines and pumps used in many industries. 
In this case, the faults are clogging and leakages. The 
SRA proposed by Gertler is further simplified and 
applied to Three-Tank System. The SRA involves two 
steps i) generation of Primary Residual Vector (PRV) 
for fault detection and ii) transformation of PRV into 
structured residual vector (SRV) for fault isolation. 

The implementation procedure of the proposed FDI 
scheme is illustrated in Figure 1. The controller in the 
system is used to maintain the process variable at its set 
point. When there is a fault in the process, its output 
differs with model output. This difference is termed as 
residual. By simply monitoring the residuals one can 
say that something is going wrong. But it is not 
possible to identify the location of the fault. So the 
residual has to be processed to enhance isolation. In 
this paper the structured residual approach is applied to 
a MIMO system to enhance fault isolation. This paper 
also considers the isolation of multiple faults that may 
occur at the same time. Theoretically, any faulty 
scenarios, e.g. one, two, or multiple simultaneous faults 
can occur. However, the probability of such scenarios 
can vary from high to low. 




Figure 1 . Block diagram representation of proposed FDI scheme 
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This paper is organized as follows. In section 2 the system 
under study which is the three-tank system is described. 
In section 3, the identification of un-measured disturbance 
variables (faults) using residual approach as reported in 
literature is explained. The proposed scheme to identify 
and estimate the unmeasured disturbance acting on the 
process is presented in section 4. In section 5 the 
simulation results are discussed. Finally the conclusions 
are drawn and scope of further work is provided in section 
6 . 

2. System descriptions 

The three-tank system considered for study [6] is shown 
in Figure2. The controlled variables are the level of the 
tankl (hi) and level of the tank3 (h3). In flow of tankl 
(fini) and in flow of tank3 (fm 3 ) are chosen as 
manipulated variables to control the level of the tankl and 
tank3. The unmeasured outflow of that is leak of tankl, 
tank2 and tank3 have been considered as fault variables 
(Li ,L 2 and L 3 ).In this paper it is assumed that the leaks 
are independent of level of the tank. 






fin 3 















^ Sl^ 








^"s3 




hi 




h2 




h3 








C } 



If 



W 



T L, ’ ^ L 3 

Figure .2. Three Tank System 
The material balance equation for the above three-tank 
system is given by 



dt SI SI ^ 

IT Vw-M) -4 

, 4 

dt S3 S3 ^ S3 V 3 



( 1 ) 

( 2 ) 

( 3 ) 



The steady state operating data of the Three-tank system 
is given in Table 1. The state space model of the three tank 
system around the operating point which are given in the 
table is 
x = Ax + Bu 

y = Cx + Du 
Where 



C = 



D- 



0 0 



0 0 
0 0 
0 0 



When fault (leak) occurs the state space model is given 
by 

x — Ax + Bu + B y f 
y = Cx + Du 

64.935 0 0 

0 64.935 0 

0 0 64.935 



B f = 



3. Fault [leak] detection using residual 
generator 

Residuals are generated from the observable variable of 
the monitored plant, that is, from the command values 
of the controlled inputs and the outputs [6]. 

Table. 1 Steady state operating data 



hi, h2,h3 in m 


0.7, 0.5, 0.3 


fin! and fm 3 in ml/sec 


100 


Outflow coefficient 

(Azl, Az2, Az3) 


2.251e-5,3.057e-5,2.307e-5 


Area of tank (Sl-S3)in 
m 2 


0.0154 


Li,L 2 ,L 3 in ml/sec 


0 


Acceleration due to 
gravity in m/sec 2 


9.81 



Ideally, the residuals should only be affected by the 
faults. However, the presence of disturbances, noise 
and modeling errors also causes the residuals to 
become nonzero and thus interferes with the detection 
of faults. Therefore the residual generator needs to be 
designed so that it is maximally unaffected by these 
nuisance inputs, which means that it is robust in the 
face of disturbance, noise and model errors. 

Structured residual are so designed that each residual 
responds to a different subset of faults and insensitive 
to the others. When a particular fault occurs, some of 
the residuals do respond and others do not. Then the 
pattern of the response set, the fault signature or fault 
code, is characteristic of the particular fault. 





0.00723996 


0.00723996 


0 


Li 


IX 2 

1 0 


3 ~^l 

0 


A = 


0.00723996 


-0.01707086 


0.0098303 


Example for fault code;- L 2 


0 


1 


0 




0 


0.0098309 


-0.012540553 


















l 3 


0 


0 


1 



B = 



64.935 
0 
0 



0 

64.935 

0 



0 
0 

64.935 



The above fault code implies that fault LI affects only 
residual R1 like L2 affects R2 and L3 affects R3.In 
order to perform detection and isolation of set of faults, 
structured residuals can be used. The so called signature 
code describes the subset of residuals which react to 
each fault. Since the levels of all three-tank are assumed 
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to be measurable, there will be three residuals 
corresponding to each of the three tanks When there is a 
fault, all the three residuals get affected. By simply 
monitoring the residuals it is possible to predict the 
change in behavior of the system from normal. But it is 
not possible to identify the location of the fault. So the 
residual has to be transformed to enhance isolation. 

Let the plant output is given by 

Y p (s) = G(s)U(s) + G f (s)L(s) (4) 

Where 

G(s) ^Transfer function under normal conditions 
G f (s) = Fault Transfer function 
U(s) =input to the system 
L(s) =fault input (leak) in the system 
Let the output of the model be given by 
Y m (s) = G(s)U (s) 

Where 



( 5 ) 



U(s) = 



L(s)-- 



fi n \ 0) 
o 

fm 3 (s) 

L\ (■*) 

L 2 (s) 
L 3 (s) 



R,(s) = [Z(s)G f -\s)][R(s)] 

Substituting R(s) in the above expression 

^(^)=[Z(5)G f - 1 (5)][G f ( 5 )I(5)] 

R t (s) = Z(s)IL(s) 

(13) 

Where I is the identity matrix. 

R t ( s ) - Z(s)L(s) (14) 

If Z(s) is the diagonal matrix then 

R lt {s) = Z n (s)L l {s) (15) 

R 2 t (s) = Z 22 (s)L 2 (s) (16) 

R 3 t (s) = Z 33 (s)L 3 (s) (17) 

It is inferred that the first element of 7^(s) that is 
Rit(s) affected only if there is a leak in the first tank. 
The second element of 7^(s) that is R 2 t(s) affected 
only if there is a leak in the second tank. Like that the 
third element of 7^(s) that is R 3 (s) affected only if 
there is a leak in the third tank. 

The transformation matrix W(s ) is 



W(s)= 



Where 



Residual R( s) is defined as difference between process 
output and model output. 

R(s) =Y p (s)-Y m (s) (6) 

Substituting the expressions for Y p (s) and 
Y m (s) in equation (6) 

R (s) =G f (s) L (s) 



W n = 



W n 


012 


W n 


W 2l 


W , 22 


023 


w 3 1 


W 32 


033 


1 


-0.0011 



(18) 



0.0154s + 0.00011 



W l3 = 0 
Wrr 



- 0.0011 



0.0154s + 0.00011 



W 22 = 1 



-0.0015 



^iO) 




G fu 


G F\2 


G F13 


F(s) 




w 23 - 

0.01545 + 0.00011 


R 2 {s) 


= 


G F2\ 


G F22 


G F23 


L 2 (s) 


(7) 


W 3 1=0 


R:M)_ 




_^F31 


G F32 


G F33 _ 


L 3 (s)_ 




jjz -0.0015 

w 31 - nniirA 



( s ) — G F i ] (s)L, (s) + G Fl 2 (s)L 2 ( s ) + G F] 3 (s)Lj ( s ) ( 8 ) 

R 2 ( s ) = Gp 2 1 G)f] (s) + Gp 22 (s)L 2 (5) + G F23 (s)L^ (5) ( 9 ) 

R 3 ( s ) = G F 3 ] (sMs) + G F 32 (s)L 2 (s) + G F 33 (s)Ln > (s) (lb) 

The above equation is valid only in the absence of plant 
model mismatch and in the absence of state and 
measurement noise. From the above equations, it is 
evident that the presence of fault will affect all the three 
residuals. It is difficult to identify the location of fault by 
monitoring the residuals. Therefore in order to enhance 
the fault isolation it is required to transform the residuals. 
The design of transformation matrix is discussed in the 
subsequent section. 

4. Design of transformation matrix 

To transform raw residual R(s) into structured form R t (s), 
multiply R(s) with weighting matrix W(s) 

R t (s) = W(s)R(s) (11) 

Weighting matrix is chosen as to contain inverse of 
Identified fault model to cancel the effects of fault 
transfer function present in residual. 

W(s) = Z(s)G f ~\s ) (12) 

Where Z(s) is user defined matrix 
Substituting W(s) in expression (11) 



0.0154s + 0.00011 



W 33 =l 

The user specified residual specification matrix Z(s) is 
chosen in the way not disturbing fault matrix even it 
can be identity matrix and it is given by 



Z(S): 



0 



"22 



0 

0 

^33. 



(19) 



Where Z n = 



1 






Z 55 = 



0.0154s + 0.0001 12 
1 

0.0154s + 0.00026 
1 

0.0154s + 0.000193 



5. Simulation results 

The proposed FDI scheme has been implemented on a 
three-tank system and its performance is observed. The 
controlled variables are the level of tankl (hi) and 
tank3 (h3). Inflow of tankl (fmi) and tank3 (fm 3 ) are 
chosen as manipulated inputs. Outflow of tankl (Li), 
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tank2 (L 2 ) and tank3 (L 3 ) are considered as leak variables. 
The synthesis method [15,16] is used for the design of PI 
controller. The PI controllers are designed so that the 
closed loop process behaves like a first order system with 
unity gain and time constant same as the open loop time 
constant. The resulting parameter for controlling the 
height of tankl using the inflow of tankl is given by Kc = 
2.54e-04 (ml/sec/m) and Ti=222 seconds and that of 
tank3 using the inflow of tank3 is given by Kc=7.69e-4 
(ml/sec/m) and Ti=200 seconds. The process is simulated 
using the non-linear first principles model, whereas the 
FDI is based on the time invariant linearized model 
(Transfer function model). 

The closed loop behavior of the process when a leak of 
magnitude 50ml/sec introduced at time t= 3000 seconds in 
tankl is shown in Figure3. The behavior of the residuals 
is shown in Figure4. From the Figure4 one can infer the 
presence of leak in tankl affect all the three residuals. By 
simply monitoring either the Process output or the 
residual it is not possible to identify the location of the 
fault. The Structured residual output is shown in Figure 5, 
6&7. From these figures one can conclude that there is 
leak only in the first tank. 

Closed loop response of the System when leak occurs in 
all the three-tanks is shown in the Figure 8 that is, a leak 
of magnitude 50ml/sec given in tankl at t=2750sec, leak 
of magnitude 50ml/sec given in tank2 at t=4000sec and 
leak of magnitude lOOml/sec given in tank3 at t=6000sec. 
It is observed that the levels of the tank are maintained 
even though the fault occurs in the process. So simply 
monitoring the process output it is not possible to detect 
the fault. The behavior of the residuals is shown in 
Figure9.With the residual one cannot find the location of 
fault. The Structured residual outputs are shown in Figure 
10, Figure 11, and Figure 12. From these figures one can 
conclude that there is leak in the all the three tanks. 

The structured residual approach is tested for modeling 
errors that is 10% deviation in time constant is 
considered. The closed loop behavior of the process when 
a leak of magnitude 50ml/sec introduced at time t= 6000 
seconds in tank3 is shown in Figure 13. The behavior of 
the residuals is shown in Figure 14. The Structured 
residual outputs are shown in Figure 15, Figure 16, and 
Figure 17. From these figures one can conclude that there 
is leak in the third tank only. 

Structured residual approach is also tested for another set 
of controller parameters. The PI controller settings are 
obtained so that the closed loop process behaves like a 
first order system with unity gain and time constant 
10%less than the open loop time constant. The closed 
loop behavior of the process when a leak of magnitude 
lOOml/sec introduced at time t= 6000 seconds in tank3 
under the new settings is shown in Figure 18. The 
behavior of the residuals is shown in Figure 19. The 
Structured residual outputs are shown in Figure20, 
21&22. From these figures one can conclude that there is 
leak in the third tank only. 

The structured residual approach is tested for slope fault. 
Leak is introduced as a ramp starting with zero value at 
t= 4000 seconds and reaching a maximum value of 100 
ml/sec at t=5000 seconds, at which it remains constant up 



to time t<=6000 seconds after that a negative ramp 
starts and it’s reaching zero at t=8000 seconds. Closed 
loop response of the System when ramp leak introduced 
is shown in the Figure23. It is observed that the levels 
of the tank are maintained even though the fault occurs 
in the process. The behavior of the residuals is shown in 
Figure24. The Structured residual outputs are shown in 
Figure25, 26&27.From these figures it is inferred that 
the proposed scheme can identify slope fault. 




Sampling instants 



Figure. 3 Closed loop response of the System when a leak of 
magnitude 50ml/sec occurs in tankl at t=3000sec 




Figure 4. Evolution of residuals when step change in 
leak of tankl of magnitude 50ml/sec introduced at 
t=3000 seconds 



a> 

T3 







Sampling instants 



Figure 5. Structured residual 1 of the system when leak of magnitude 
50ml/sec is given in Tankl at t=3000 seconds 
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Sampling instants 



Figure 6. Structured residual2 of the system when leak of magnitude 
50ml/sec is given in Tankl at t=3000 seconds 




Sampling instants 



Figure 7. Structured residuaB of the system when leak of magnitude 
50ml/sec is given in Tankl at t=3000 seconds 




Figure 8. Closed loop response of the System when a leak of magnitude 
50ml/sec given in tankl at t=2750sec, leak of magnitude 50ml/sec given 
in tank2 at t=4000sec and leak of magnitude lOOml/sec given in tank3 at 
t=6000sec 




Figure 9. Evolution of residuals when a leak of magnitude 50ml/sec 
given in tankl at t=2750sec, leak of magnitude 50ml/sec given in tank2 
at t=4000sec and leak of magnitude lOOml/sec given in tank3 at 
t=6000sec 




Samling instants 



Figure 10. Structured residual 1 of the system when a leak of magnitude 
50ml/sec given in tankl at t=2750sec, leak of magnitude 50ml/sec given 
in tank2 at t=4000sec and leak of magnitude lOOml/sec given in tank3 at 
t=6000sec 




Figure 1 1 . Structured residual2 of the system when a leak of magnitude 
50ml/sec given in tankl at t=2750sec, leak of magnitude 50ml/sec given 
in tank2 at t=4000sec and leak of magnitude lOOml/sec given in tank3 at 
t=6000sec 




5 
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Figure 12 Structured residuaB of the system when a leak of magnitude 
50ml/sec given in tankl at t=2750sec, leak of magnitude 50ml/sec given 
in tank2 at t=4000sec and leak of magnitude lOOml/sec given in tank3 at 
t=6000sec 




Sampling instants 



Figure 16. Structured residual2 of the system when leak of magnitude 
lOOml/sec is given in Tank3 at t=6000 seconds 




Figure. 13 Closed loop response of the System when a leak of magnitude 
lOOml/sec occurs in tank3 at t=6000sec 




Figure 17. Structured residuaB of the system when leak of magnitude 
lOOml/sec is given in Tank3 at t=6000 seconds 




Sampling instants 



Figure 14. Evolution of residuals when step change in leak 
of tank3 of magnitude lOOml/sec introduced at t=6000 
seconds 




Sampling instants 



Figure 15. Structured residual 1 of the system when leak of magnitude 
lOOml/sec is given in Tank3 at t=6000 seconds 




Figure. 18 Closed loop response of the System when a leak of magnitude 
lOOml/sec occurs in tank2 at t=4000sec 




Figure 19. Evolution of residuals when step change in leak of tank2 of 
magnitude lOOml/sec introduced at t=4000 seconds 




6 
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Sampling instants 



Figure 20. Structured residual 1 of the system when leak of magnitude 
lOOml/sec is given in Tank2 at t=4000 seconds 




Figure24. Evolution of residuals when a leak of varying magnitude 
(ramp) lOOml/sec occurs in tankl at t=4000sec to t=8000 sec 




Figure 2 1 .Structured residual2 of the system when leak of magnitude 
lOOml/sec is given in Tank2 at t=4000 seconds 




Sampling instants 



Figure22. Structured residuaB of the system when leak of magnitude 
lOOml/sec is given in Tank2 at t=4000 seconds 




sampling instants 



Figure.23 Closed loop response of the System when a leak of varying 
magnitude (ramp) lOOml/sec occurs in tankl at t=4000sec to t=8000 sec 




Figure 25. Structured residual 1 of the system when leak of varying 
magnitude (ramp) lOOml/sec occurs in tankl at t=4000sec to t=8000 sec 




Sampling instants 



Figure 26. Structured residual 1 of the system when leak of varying 
magnitude (ramp) lOOml/sec occurs in tankl at t=4000sec to t=8000 sec 




Figure 27. Structured residual 1 of the system when leak of varying 
magnitude (ramp) lOOml/sec occurs in tankl at t=4000sec to t=8000 sec 



7 
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6. Conclusions 

The performance of the proposed scheme has been 
evaluated on a three- tank process for leak in the tanks. 
The proposed FDI scheme can provide fault information 
even when there is simultaneous change in more than one 
leak. It should be noted that the proposed method is 
independent of the controller design. From the structured 
residuals, the magnitude of leak and time of occurrence of 
leak are also found. And one can conclude that the 
estimated magnitude and time of occurrence of the leak 
variable (fault) are close to the true value. So from the 
proposed method one can identify the fault as soon as it 
occurs in the process. The proposed method is found to be 
robust to plant model mismatch. 
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Abstract 

In this paper, a transient model of the faulty machine is 
developed. The model is referred to a three phase stator 
winding, while the rotor has been represented by all the 
meshes allowing for the representation of various faults. 
The model is based on coupled magnetic circuit theory 
and incorporates non sinusoidal air-gap magneto motive 
force (MMF) produced by both stator and rotor, therefore 
it will include all the space harmonics in the machine. 
Simulation and experimental results were first used to 
identify low and high frequency spectral components 
created by rotor anomalies in the stator current spectrum. 
This paper then presents an approach for motor rotor 
fault diagnosis using neural networks. The proposed 
methodology has been experimentally tested on a 
5.5Kw/3000rpm induction motor. The obtained results 
provide a satisfactory level of accuracy. 

Keywords : Induction machines, coupled magnetic 

circuit, rotor failures, Neural networks, diagnostic. 

1. Introduction 

The use of induction motors in today’s industry is 
extensive, and the motor can be exposed to different 
hostile environments, misoperations, and manufacturing 
defects. Internal motor faults (e.g., short circuit of motor 
leads, intertum short circuits, ground faults, bearing and 
gearbox failures, broken rotor bar and cracked rotor end- 
rings, as well as external motor faults (e.g., phase failure, 
asymmetry of main supply and mechanical overload), are 
expected to happen sooner or later [1]. Furthermore, the 
wide variety of environments and conditions that the 
motors are exposed to can age the motor and make it 
subject to incipient faults. These incipient faults can lead 
to motor failure if left undetected. 

With proper system monitoring and fault detection 
schemes, the costs of maintaining the motors can be 
greatly reduced and the availability of these machines can 
be significantly improved. 

Many researchers have focused their attention on 
incipient fault detection and preventive maintenance in 
recent years. There are invasive and noninvasive methods 



for machine fault detection, the noninvasive methods are 
more preferable than the invasive methods because they 
are based on easily accessible and inexpensive 
measurements to diagnose the machine conditions 
without disintegrating the machine structure. 

Recently, artificial intelligence (AI) techniques have been 
proposed for the noninvasive machine fault detection 
[2,3]. These Al-based techniques include expert systems, 
neural network, fuzzy logic and combined techniques. 
These Al-based techniques have numerous advantages 
over conventional diagnostic approaches. In general, they 
can lead to improved performance when properly tuned, 
they are easy to extend and modify, and they can be 
easily made adaptive by the incorporation of new data or 
information, as they become available. 

In the Al-based systems, several quantities are utilized as 
input signals: stator currents and voltages, magnetic 
fields and frame vibrations, etc. In general, stator currents 
and voltages are preferred because they allow for the 
realization of non-invasive diagnostic systems 
The diagnostic procedure, based on the current signal, 
can be organized in the following steps [2]: 

■ Choice of the failure to be considered, 

■ Definition of cause-effect relationships, 

■ Computation of the diagnostic indexes linked to the 
fault extent. 

Filippetti[4] note that the AI techniques to induction 
machines diagnostics should correspond to wide 
industrial applications, but plant operators look at these 
tools with a certain degree of suspicion, due to their 
complexity and their cost. To implement a complete 
diagnostic system for industrial applications, it is 
fundamental to: 

■ Employ a minimum configuration intelligence, 

■ Develop a knowledge base transferable from 
machine to machine. 

With regard to the first point, in connection to the 
structure of a diagnostic system, the best idea appears to 
be a global neural solution. 

In order to develop technologies for motor fault detection 
and diagnosis, it is important to show the proposed 
theory, schemes, feasibility, and limitations. In the 
research stage, we require a controllable environment so 
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that we know what types of faults are induced and what 
types of motor performance are resulted. 

The objective of this paper is to develop a model which is 
capable to predict the performance of induction machines 
under broken rotor bars and end-ring faults then a 
supervised neural network is applied to motor rotor faults 
detection and severity evaluation. Neural network is 
trained and tested using measurement data of stator 
current spectra. 

The remainder of the paper is organized as follows: 
Section (2) focuses on the complete model of a squirrel 
cage induction machine. In section (3) we discuss the 
modeling of rotor failures. In Section (4) simulation 
results are present. In Section (5) we describe the 
experimental setup and results. In section (6) the neural 
networks-based detection system is implemented. We 
present in section (7) results using neural network. 
Section (8) concludes this paper. 

2. Induction machine modeling 

This model follows the coupled magnetic approach by 
treating the current in each rotor bar as an independent 
variable. The effect of non-sinusoidal air-gap MMF 
produced by both the stator and the rotor currents have 
been incorporated into the model. This is done by use of 
the winding function approach. 

Our analysis is based on the following assumptions [5]: 
Symmetric machine, uniform air-gap, negligible 
saturation and insulated rotor bar. 

The stator comprises of three phase concentric winding. 
Each of these windings is treated as a separate coils. The 
cage rotor consists of n bars can be described as n 
identical and equally spaced rotor loop [6,7]. Each loop is 
formed by two adjacent rotor bars and the connecting 
portions of the end-rings between them. Hence, the rotor 
circuit has n+1 independent currents as variables. The n 
rotor loop currents are coupled to each other and to the 
stator windings through mutual inductances. The end- 
ring loop does not couple with the stator windings [7], it 
however couples the rotor currents only through the end 
leakage inductance and the end-ring resistance. 

2.1. Stator voltage equations 

The stator equations for the induction machine can be 
written in vector matrix form as: 
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dt 
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where 
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The matrix [R s ] is a diagonal 3 by 3 matrix which 
consists of resistances of each coil. 

Due to conservation of energy, the matrix [Z 5 J is a 
symmetric 3 by 3 matrix. The mutual inductance [ L sr \ 
matrix is an 3 by n matrix comprised of the mutual 
inductances between the stator coils and the rotor loops. 
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Where L srij is the mutual inductance between the stator 
phase i (i=l,2 or 3) and the rotor loop j ( j=l,2 ... n) and 
L srie the mutual inductance between the stator phase i 
(i=l,2 or 3) and the end-ring. 

2.2. Rotor voltage equations 

Given the structural symmetry of the rotor, it is 
convenient to model the cage as identical magnetically 
coupled circuits. For simplicity, we assume that each 
loop is defined by two adjacent rotor bars and the 
connecting portions of the end-rings between them [7]. 
For the purpose of analysis, each rotor bar and segment 
of end-ring is substituted by an equivalent circuit 
representing the resistive and inductive nature of the 
cage. Such an equivalent circuit is shown in Figure 1. 
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Figure 1. Rotor cage equivalent circuit showing rotor loop currents and 
circulating end ring current. 



From Figure 1, the voltage equations for the rotor loops 
can be written in vector matrix form as: 
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where 
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ring leakage inductance and L ki is the mutual inductance 
between two rotor loop. 

2.3. Calculation of torque 

The mechanical equation of the machine is [5]: 



/ 

r 

In case of a cage rotor, the rotor end ring voltage is 
and the rotor loop voltages are V rk =0, k = 1,2. . .n. 
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Ke=0, 



The loop equation for k th rotor circuit is: 



C — C —J 



em r 



dQ. 

m_ 

dt 



(14) 



Where C em is the electromagnetic torque, C r is the load 
torque , J is the inertia of the rotor and f2 m is the 
mechanical speed. 
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The voltage equation for the end-ring is: 



V =0=-RI ,-RI 

re e r i e rl 



d 0> 

• — R I -\-tiR I H — — 

e rn e re 



where 6 is the angular position of the rotor and p denotes 
the number of motor pole pairs. 

(11) The electromagnetic torque is given by the following 
equation [8]: 



Where R b is the rotor bar resistance, R e is the end-ring 
segment resistance and O rk is the total flux linked by the 
Ml loop. 

Since each loop is assumed to be identical, the equation 
(10) is valid for every loop. Therefore the resistance 
matrix [R r ] is a symmetric (n+1) by (n+1) matrix given 
by: 
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( 12 ) 

Due to the structural symmetry of the rotor, [ L rr ] can be 
written in matrix form (13). 

[4 
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(13) 

where L kk is the self inductance of the Ml rotor loop, L b 
is the rotor bar leakage inductance, L e is the rotor end- 



C =PI 




(16) 



2.4. Calculation of inductances 

It is apparent that the calculation of all the machine 
inductances as defined by the inductances matrices in the 
previous section is the key to the successful simulation of 
an induction machine. 

The model must take into account the geometric 
construction of the machine and then will include all the 
space harmonic. 

These machine inductances are conveniently calculated 
by means of winding functions. This method assumes no 
symmetry in the placement of any motor coil in the slots. 
According to winding function theory, the mutual 
inductance between two windings i and j in any electric 
machine can be computed by the following 
equation [8,9,10]: 

L ( 6 ) = toll 2 \ Si (6,<p)s/6,<p)d<p ( 17 ) 

go 7 

where ju 0 = 471. 1 0" 7 H/m, g is the air gap, 6 is the angular 
rotor position, r is the average radius of the air gap, / is 
the active length of the motor, ^?is a particular point 
along the air gap and £i(0,cp) is called the winding 
function and represents the magneto motive force (MMF) 
distribution along the air gap for a unit current flowing in 
winding i. 

For simulation purposes a three phase, 4Kw, 50hz, 4 
pole, 3 80/220 V squirrel-cage induction motor will be 
treated in this paper. The machine has a stator comprises 
of three phase concentric winding with 36 slots, 28 rotor 
bars and two coils per phase. Each coils is constituted of 
three sections having N s turns in series (N s =26 for the 
machine under test). Figure 2 shows the MMF 
distribution produced by 1A of current through the stator 
phase " 1 " obtained by summing all the winding function 
of the stator phase coils [10]. Note that the both stator 




11 




Automatic Control and System Engineering Journal, ISSN: 1687-4811, Volume 7, Issue 2, ICGST, Delaware, USA, Nov. 2007 



phases "2" and "3" produce a similar MMF distribution 

but shifted by — and — . 

3 3 

The MMF distribution produced by 1A of current 
through a rotor loop can only take two values depending 
on whether we are inside or outside the loop. The angle 

between two adjacent rotor bars is a = —, the MMF 

14 

distribution produced by 1 A of current through the first 
rotor loop, is shown in Figure 3. 

The mutual inductance between stator and rotor branches 
will be a function of the rotor position angle, 6. Figure 4 
gives the mutual inductance (L srll (6)) between the stator 
phase "1" and the rotor loop "1". Note that the mutual 
inductance between the phase "2" and the rotor loop "V 
is L srll (0) but shifted to the right by 6y where y is the 
angle between two stator slots. Mutual inductance 
between the phase "1" and the rotor loop "2" is L srll (6) 
but shifted to the left by a, where a is the angle between 
two rotor bars. 




Figure 2. MMF distribution of the stator phase "l" 



a * 










2n 




a 


n 


2n j 


a . 






2ti 











Figure 3. MMF distribution of the rotor loop " 1 " 
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bar or end-ring faults. In the case of m fully broken bars 
or end-ring segments, m loop equations are removed. 

4. Simulation results 

To validate the proposed model, a functional schema of 
the induction machine was developed on the MATLAB- 
SIMULINK platform. 

The machine is first simulated under healthy condition. 
Figure 5 shows the instantaneous electromagnetic torque, 
speed and the phase current of the machine during a start 
up with a balanced sinusoidal voltage supply followed by 
the application of a load (C r = 1 INm) at instant t=0.8s. 
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Figure 4. Mutual inductance between the stator phase " 1 " and rotor 

loop "1" 



Figure 5. Torque, Speed and stator current in phase "l"(top to 
bottom). Normal machine 



3. Modeling rotor bars faults 

Broken bar faults can be incorporated in the healthy 
machine model by assuming that bars or end-ring 
segments are fully broken [6]. Thus, one can remove the 
loop equations corresponding to the broken bars or the 
end-ring segments from (7) in order to model the broken- 



Next the rotor failures were simulated. Figure 6 Shows 
the instantaneous electromagnetic torque, speed and the 
phase current of the machine with four broken rotor bars 
and one end-ring, during a start up with a balanced 
sinusoidal voltage supply followed by the application of a 
load (C r = 1 INm) at instant t=0.8s. We can easly see that 
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the effect of four broken rotor bars and one end-ring is 
very important, (e.g., the acceleration time under rotor 
asymmetry is larger than under healthy machine). 




Time (s) 




Figure 6. Torque, speed and stator current in phase"l"(top to bottom). 
Machine with four broken rotor bars and one end-ring 

4.1. Analysis of steady state operation 

The stator current signal at steady state for the loaded 
machine is transformed by the Fast Fourier Transform 
(FFT) into signal in the frequency domain to generate the 
power spectral density (PSD). The spectrum generated by 
this transformation includes only the magnitude 
information about each frequency component, which can 
be analyzed and processed easier than signal in the time 
domain. 

4.1.1. Spectrum analysis in the bandwidth 
[25Hz-75Hz] 

Figure 7 reports stator current spectrum around 
fundamental for a normal machine and a machine with 
four broken bars and one broken end-ring. We can see 



that rotor anomalies induce some harmonic components, 
given by [4,1 1]: 

f bc = (\±2ks)f s ,£=1,2,3, (18) 

where are frequencies associated with the broken bar, 
s is the slip and f s is the supply frequency. 




frequency (Hz) 




Figure 7. Simulated stator current spectrum for healthy machine and 
for the case of four broken rotor bars and one end-ring ( top to bottom) 

4.1.2. Spectrum analysis in the bandwidth 
[lOOHz-lOOOHz] 

There are other spectral components that can be observed 
in the stator line current due to broken rotor bar fault. The 
equation describing these frequency components is given 
by Nandy [11] 




\ 



J 




f S 



(19) 



± 

where, f hbc are detectable broken bar frequencies; 

1 = 3 , 5 , 7 , 9 , 11 , 13 , 



Figure 8 shows the simulated plot of the stator currents 
spectrum around 5th time harmonic with rotor failures, 
we can visualize the presence of the principal 

+ 

components f hbc . We notice the presence of additional 

frequency components which are spaced by 2sf s . This can 
be verified on space harmonic 7, 11, 13 ... 
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Then, it is necessary to complete equation (19) to take 
into consideration these harmonics because they are 
indicative of failure presence in the rotor cage. The new 
equation will be: 



fh 



hbc 



vL 



(l-.sj± (] +2/?),s 



fs 



( 20 ) 



where, f hbc is detectable broken bar frequencies; 
-^=3, 5, 7, 9,1 1,13, And /?= 0,1, 2, 3, 




frequency (Hz) 



Figure 8. Simulated stator current spectrum for faulty machine around 
the 5th time harmonic 

5. Experimental setup and results 

The characteristics of the 3 -phase induction motor used 
in our experiment are listed in Table 1. The needed load 
of the induction motor was established by connecting the 
test motor to an eddy current brake via a flexible 
coupling (Figure 9). In order to allow tests to be 
performed at different load levels, the brake DC supply 
current is controllable. 




Figure 9. View of the experimental setup 



A current Hall effect sensor (Honeywell type CSNA111) 
was placed in one of the line current cables. The stator 
current was sampled with a 4 KHz rate and interfaced to 
a pentium PC by an ARCOM(APCI-ADADIO) 
acquisition board. The APCI-ADADIO is a 32-bit PCI 
local bus board, which provides 8 differential or 16 
single-ended multiplexed 12-bit ADC channels. 



A MATLAB Application Program Interface (API) is 
built to allow The ARCOM A/D PC card to communicate 
directly with MATLAB. 

The motor was tested with the healthy rotor and a faulty 
rotor with two broken bars. The bars were broken by 
drilling holes through them. 



Description 


Value 


Power 


5.5 kW 


Input Voltage 


220/380 V 


Full load current 


20.6/1 1.9A 


Supply frequency 


50 Hz 


Number of poles 


2 


Number of rotor slots 


28 


Number of stator slots 


36 


Full load speed 


2875 rpm 



Table 1. Induction motor Characteristics used in the experiment 



As predicted by simulation, Figures 10, 11 and 12 show 
the experimental results with healthy and faulty machine 
around fundamental. We can see the amplitude of the 
sidebands components f bc according to equation (18) in 
the faulty machine current increase over their counterpart 
in the case of the healthy machine. 

Experimental results show significant changes around the 
5th time harmonic (Figure 13). According to equation 
(20), space harmonic spaced by 2sf s can be clearly seen. 




frequency (Hz) 

Figure 10. Experimental stator current spectrum around fundamental of 
the healthy machine 




frequency (Hz) 



Figure 11. Experimental stator current spectrum around fundamental of 
the faulty machine with one broken bars 
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Figure 14. : Block diagram of the experimental NN-based detection system 




Figure 12 Experimental stator current spectrum around fundamental of 
the faulty machine with two broken bars 




Figure 13. Experimental stator current spectrum around 5th time 
harmonic of the faulty machine with two broken bars 



6. Description of the experimental Neural 
Networks-based detection system 

Stator current analysis is based on the principle, as the 
fault progresses, its characteristic spectral components 
continue to increase over time. Therefore, the correlation 
between the amplitude of these components and the fault 
extent is the issue of the diagnostic system. 

In this section, real motor stator current data were 
collected to verify the feasibility of applying a neural 
network to diagnose rotor fault. 

Figure 14 describes the basic fault detection system in the 
frequency domain. First stator current signal are collected 
from Hall effect sensor at an equal sampling time. Then 
the signals are transformed by the Fast Fourier transform 
(FFT) into signals in the frequency domain, which can 
provide salient features for the diagnosis of rotor 
condition. The spectrum generated by this transformation 
includes the magnitude information about each frequency 
component. Signal noise that exists in the spectrum is 
reduced by averaging and windowing. 

The ability to average a series of measurements is useful 
to discriminate between noise and components that are 
actually part of the signal. This ensemble-averaging 
technique is very effective for determining the frequency 
content of a signal buried in a random noise. 

Finally, the rotor faults frequency signature is extracted 
as inputs of neural networks. Through supervised training 
with inputs and outputs, the learned neural networks can 
detect rotor faults and also the extent of the faults. 

6.1. Artificial Neural Networks (ANN) 

Artificial neural networks (ANN) are used to recognize 
and classify complex fault patterns without much 
knowledge about the system they deal with. They are 
designed to mimic the human nervous system using 
massively parallel nets composed of many computational 
elements connected by links with variable weights. Of all 
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the ANNs, the multi-layer feedforward artificial neural 
network, trained using the back-propagation algorithm, is 
the most commonly and flexibly used [12]. Its typical 
architecture, which contains the input layer, a number of 
hidden layers and the output layer, is shown in Figure 15. 
Initially, an input pattern is applied to the input layer, 
then signals propagate forward, level by level, through 
the hidden layers and the output layer produces the 
network output. Each neuron in a layer (except the ones 
in the input layer) sums up their input value (x) with 
interconnection to its associated weight (w) and 
corresponding bias ( b ) to form its own scalar output 
{Net). The bias is added to shift the sum relative to the 
origin. The Net passes through a non-linear activation 
function (/), providing an output value (Y). 

The activity level of the jth neuron in layer / is obtained 
as [13]: 






i 



( 21 ) 



where Oj is the activity level (output) of the yth neuron, 
Netj is the input of the yth neuron,^ is the activation 
function of the yth neuron, w Jt is the connection weight 
from the ith neuron to the jth neuron, v 7 is the activity 
level of the zth neuron in the prior layer and bj is the bias 
term of the yth neuron. 



Input layer Hidden layers Output layer 




= * k -y t ) (23) 

k= 1 v y 

Where 4 represents the desired output of the Ml neuron 
in the output layer. 

The incremental change of weight from the zth neuron to 
the yth is computed by [14]: 

A w (t+1 ) = rj 8 . 0 . + a A w (t) (24) 



s k=( t k-°M Net k) 



(25) 




(26) 



where A Wj t {t) is the incremental change in the weight wy 
at time t , t k is the desired output of the kth neuron in the 
output layer, rj is the learning rate and a is the 
momentum. 

The momentum term ( aAw..(t ) ), included in the weight 

update equation to try to avoid a local minimum [1]. 
Equation (25) holds for the Ml neuron in the output layer, 
and (26) holds for the mth neuron in the hidden layer. 



6.2. Data normalization 

In order to improve the neural network performance, the 
data must be well-processed and properly-scaled before 
inputting them to ANN[13]. The normalization process 
was required to restrict the range of the patterns for input 
into the neural networks. In this study, the inputs to the 
neural network are normalized between [-1,+1] by 

equation (27): 

r \ 



x = 



x-x 

m 

x-x 



xl,8-0,9 



(27) 



Figure 15 . Basic structure of a multi-layer feedforward artificial 
neural network 

The active function defines the output of a neuron in 
terms of the activity level at its level. There are three 
basic types of activation functions: threshold, piece-linear 
and sigmoid [3]. The sigmoid function is by far the most 
common activation function. In this paper the tan- 
sigmoid transfer function is used, which generates 
outputs between -1 and +1. Its definition is expressed in 
equation (22): 



Where x max is the maximum magnitude of the input 
vector andx min is the minimum magnitude of the input 
vector 

6.3. Machine slip computation 

The current spectrum components depend on the machine 
speed or slip. Therefore to avoid the need to measure the 
machine speed using conventional transducers, we 
choose to calculate the rotor slot harmonics given by the 
following equation [4] : 



Wet)=—^mr W (22) 

l + e l + e 

The backpropagation algorithm is used to train the 
network [1]. The connection weights are iteratively 
adjusted so that the error between the network output and 
the desired output (target) for a given reference input is 
minimized. The error goal is expressed as the sum of the 
squared error (SSE), calculated as follows: 




v y 



(28) 



where P is the pole pair number, f s is the supply 
frequency and N r is the rotor slot number. 

These components are always present in the stator current 
spectrum with healthy and faulty rotor. Figure 16 shows 
these components for full load condition of the motor 
with two broken rotor bar fault. However, the analyze of 
this spectrum reports that the detection of the first slot 
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harmonic line at frequency^. is easy and represents the 

component with the maximum amplitude in this 
bandwidth. 




frequency (Hz) 



Figure 16. Experimental stator current spectrum in the frequency 
bandwidth [1200 Hz - 1450 Hz] 

7. Results using the neural network 

Experimental data were collected for each operating 
condition of the motor (healthy, rotor with one and two 
broken bar) under four different load conditions. The 
load conditions of the motor are 25%, 50%, 75% and full 
load respectively. These load condition percentages are 
determined according to the motor nameplate information 
given in Table 1. A total set of 180 sets of current was 
collected, and two thirds of them were used to train the 
neural networks. The rest of the data set was used to test 
the network’s performance. 

The Neural Network Toolbox used of MATLAB was 
employed to train the artificial neural networks for this 
investigation. A three layered feed-forward neural 
network with backpropagation algorithm was used to 
perform the desired analysis. The network topology is as 
follows: 

■ The NN input layer is constituted by a set of 14 neurons 
in order to consider 10 spectrum components centred 
around the fundamental, 2 spectrum components around 
5th harmonic and 2 spectrum components around 7th 
harmonic. We have considered the fundamental harmonic 

k 

given by equation (20) when r|=0 and — respectively 
equal to 5 and 7. 

■ A proper selection on the number of hidden neurons has 
significant effect on the network performance. In order to 
demonstrate the performance of the FFNN, the number of 
hidden neurons is varied to find the optimal design. The 
different numbers of hidden neurons applied in the 
verification are 5, 6, 7, 8, 9, 10,11 and 12. 

■The output layer of the neural network has three 
neurons representing different states of the machine. 
Since the FFNN type neural network belongs to 
supervised learning, it needs a teacher to lead it in order 
to achieve the determined goal. The expect target vectors 
were defined as three different pattern: Tl=[l -1 -1] T , 
T 2 =[-l 1 -1] T , T 3 =[-l -1 1] T for healthy, faulty with one n 
be extended broken bar and faulty with two broken bar 



respectively. 1 is set as the correct class and -1 for all the 
other classes. 

After successful training, the network was used to 
distinguish between the rotor states. The test data was 
unseen by the neural network. A learning rate of 0.01 and 
momentum of 0.95 was selected for all cases based on 
the experience. Optimum results were obtained with 1 1 
hidden neurons. The accuracy is 100% for the training 
data, and 96.66% for the testing data, the performance of 
the FFNN is summarized in Figure 17. 




Figure 17. Accuracy of training and testing data at different number of 
hidden neurons 

8. Conclusion 

A detailed model of a squirrel-cage induction machine 
has been developed. In order to simulate cage defects, the 
machine was modelled as a group of coupled magnetic 
circuits by considering the current in each rotor bar as an 
independent variable. The model can simulate the 
performance of induction machines during transient as 
well as at steady state, including the effect of rotor faults. 
Simulated and experimental data were first used to study 
rotor faults cause-effect relationships in the stator current 
and the current spectrum signature. Experimental stator 
current data were then applied to verify the feasibility of 
applying a feed-forward neural network with 
backpropagation algorithm to diagnosis rotor faults. The 
results show that FFNN is not only able to automatically 
identify rotor status, the neural network is also designed 
to determine the extent of faults. This is very useful for 
early fault detection and preventive maintenance. The 
proposed methodology has been tested on a 
5.5Kw/3000rpm squirrel-cage induction motor. 
Experimental tests have led to results with a satisfactory 
level of accuracy greater than 96%. 

Future work consists to test this Al-diagnosis rotor faults 
when the motor is connected to a motor speed drive. A 
Self-Organizing Map (SOM) could be used for initial 
classification of motor faults, as an introductory step. 
Then, the diagnosis problems of induction motors will be 
extended in the case of rotor and rolling bearing faults. 
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Abstract 

By the use of elegant operational properties of Fourier 
series, a direct computational algorithm for evaluation the 
optimal control and trajectory of linear singular systems 
with a quadratic cost functional is developed in this 
paper. The state variable, state- rate, and the control 
vector are expanded in Fourier orthogonal basis with 
unknown coefficients. The technique is straightforward 
and very convenient for computer programming. An 
illustrative example is given to demonstrate the 
applicability of the proposed method. 

Keywords: Optimal control; singular systems; Fourier 
series; operational matrix 

1. Introduction 

Various orthogonal functions have been used to solve the 
problem of optimal control in recent years. Special 
attention has been given to applications of Walsh 
functions [1], block-pulse functions [2], Laguerre 
polynomials [3], Chebyshev polynomials [4], Legendre 
polynomials [5], and Fourier series [6]. 

The main characteristic of this technique is that it reduces 
these problems to those of solving a system of algebraic 
equations; thus greatly simplifying the problem and 
making it computationally plausible. The approach is 
based on converting the underlying differential equations 
into integral equations through integration, 
approximating various signals involved in the equations 
by truncated orthogonal series, and using the operational 
matrix of integration to eliminate the integral operations. 
Singular systems, also commonly called generalized or 
descriptor systems in the literature appear in many 
practical situations including engineering systems, 
economic systems, network analysis, and biological 
systems (see e.g. [7, 8, 9]). In fact, many systems in the 
real life are singular in nature. They are usually 
simplified as or approximated by non- singular models 
because there is still lacking of efficient tools to tackle 
problems related to such systems. Significant results have 
been reported in diverse respects such as solvability 
[10,11], freqency-domain design [12], multiple-integaral 



observers [13] and integarted fault estimation and fault 
tolerant design [14]. Fourier series have been widely used 
to study the optimal control of linear systems with 
quadratic performance index. Campbell [15], Lewis and 
Mertzios [16], Marszalek [17], Paraskevopoulos [18], 
and Trzaska [19], applied the orthogonal functions for 
analysis of linear singular systems. Palanisamy and 
Balachandran [20, 21], and Balachandran and 

Murugresan [22], applied the single-term Walsh series 
(STWS) method to the analysis of linear singular system. 
This technique is also used to analyse the electronic 
circuits [23]. In this paper we extend the Fourier series 
technique to study the problem of optimal control of 
linear singular systems with quadratic performance. The 
state variable x(^) , state rate x{t} and control 

variable u(t) are expanded by Fourier series with 
unknown coefficients. 

Also the properties of Fourier series are used to relate the 
unknown coefficients of x(^) to the coefficients ofx(^) . 
Using this method the dynamics of the system are 
converted to a system of algebraic equations. The 
necessary conditions of optimality are imposed to solve 
the quadratic programming by using the Lagrange 
multipliers. 

The reminder of the paper is organized as follows: 
section (2) is about the preliminaries and problem 
statement; section (3) focuses on the main result, section 
(4) demonstrate a benchmark example. Finally section (5) 
is the conclusion. 

2. Preliminaries and Problem Statement 

2.1 Properties of the Fourier series 

Suppose a signal g(t) has finite energy over the 
interval [0, L] . Then g(t ) can be represented by a Fourier 
series as follows 

/ N ^ / 2n7it . * . , InTit. 

g(t) = go +2^Sn c °s(— — ) + g„ s m ( — - — ) (1) 

n=l ^ -Z-/ 
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* 

Where the Fourier coefficients g n , g n are given by 
1 L 

g 0 = — j g ( t)dt (2a) 

L 0 

2 p 2n7rt 7 . . _ 

= — fg(0cos — — dt, n = 1,2,3,... (2b) 

Z J 0 £ 

* 2 r , x . 2n7rt 7 1 _ . 

g n = ~T Jg(O sm — - — « = 1,2,3,... (2c) 

Z o z 

We assume that the derivative ofg"(^)in equation (1) is 
described by 

m = m =fo+ z{/,pos(^) + ./; sin^j (3) 

An approximation to /*(£)andg(£) are obtained by 
keeping those terms in (1) and (3) for which 1 < n < r and 
discarding others. This may be written in a compact form 
as 



git) = gA it) + £ \g,A it) + gJ (0) = g t ®(0 (4) 

n= 1 

fit) = fA it) + £ {f A it) + /„>* (0} = (5) 

n = 1 

Where 

<?=[& a a - a a a - sT 
F=[f c f, A - /, / £ - XT (« 

®(0=M A A - A A A - if 



With 

/ 2 72 7?/ A 1 o o 

(p n = cos — - — , n = 0,1, 2, 3, . . . 

/* • 2fl TCt i o ^5 

(p n = sin— — , 27 = 1,2,3,... 

The elements of ^( 7 ) are orthogonal in the 



L 

Interval te(0,L). i.e. J (/) n (/) m {t)dt 
0 



0 77 ^ 777 

< L 

— n = m 

{ 2 



The integration of vector O(^) defined in equation (6) 
can be approximation by 



Z=Z 



- 0 0 ... 0 0 



1 -1 
n In 

0 0 0 ... 0 0 — 0 

In 



-1 -1 



0 0 0 
0 0 0 

— — 0 

2 n 2n 

— 0 — 

4 n 4 n 



— 00 
.2 rn 



0 0 0 0 
0 0 0 0 
0 0 0 0 
0 0 0 0 



0—00 
2 rn 



(r-\)n rn 
0 0 



1 



%r-X)n 

0 



0 

J_ 

2r;r 

0 

0 



0 0 



( 8 ) 



The cross-product 


of two 


Fourier series vectors is 




A 2 


<t>A 


... <t>A 


<t>A * • 


.. m: 




<h<k 




... <j>A 


Ml ■ 


- m: 


mfit)= 


<t>A 


<t>A 


- $ 


<t>A • 


- ¥, 




tfA 


<fA 


... (/\(/> r 


/ • 






/A 


<fr<t 1 


- 


U ■ 


•• <t>£ 



(9) 




By using 

\jj m it)dt =\ L o f n f m {t)dt = 

J 0 4>l C t)dt = L 



Where m and n can assume any of the values 1,2,3,. . ..we 
can easily obtain 



1 



0 



D 



= ^A'(t)dt = L 



( 10 ) 



£<d (t')dt' = p<j> (o 



(7) 



Where P is an (2r + 1) x (2r + 1) operational matrix of 
integration and is defined as [24]. 



Where D is a (2r + 1) x (2r + 1) matrix. Thus D is a 
diagonal matrix and the simple structure of this matrix 
plays important role in the direct method for solving 
optimal control problem. Further, using equation (3) we 
get 



git)-g(0) = \‘ 0 fit')dt' 
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And from equation (4), (5) and (7), using g(0) = K T Q(t), 
we have 

(G t -K T )0(t) = F T P0(t) 

Where K T = [g(0) 0, ... 0, 0* ... O'Jthus we get 

(G t -K J ) = F t P (11) 

2.2 Problem statement 

Consider linear time invariant (LTI) singular system 



Ex{t) = Ax(t) + Bu{t ) 

x(0) =x 0 



( 12 ) 

(13) 



Where x e R p is the state of the system, u <= R z is the 
input vector of the system, i? e R pxp is assumed to be 
singular with 0 < rank ( E ) < p ,and A,B are real 
constant matrices of appropriate size. The problem is to 
find the optimal control u(t ) and the corresponding state 

trajectory x(t\ 0 <t<L, satisfying equation (12)- 
(13) while minimizing the cost function 



^x\L)Sx(L) + ^^ ]dt ll4) 



J = 

2 



In equation (14), T denotes transposition, S,Q and R 
are matrices of appropriate dimensions, S and Q are 
symmetric positive semi-definite matrices and R is a 
symmetric positive definite matrix. 

3. Main result 

3.1. Approximation of the singular system 

By expanding X p (t\ u z (t) and x p (t ) in Fourier 
series of order r , we determine the following 
approximated values, i.e. for P = 1,2,3,... p, and 

Z = 1,2,3,... z we get 



x p (0 = Vo (0 + £ [dp^ it) + a Pi 0* ( t ) } = A t O(/) (15) 

(=1 

U Z (0 = ^Z0^0 (0 "I" ^ fasti ^ ^Zifii (0 } = ® ®(0 (16) 

i=l 

i P (0 = c^o(0+ZK^(0+44*(0 } = c T o(0 ( 17 ) 

/=l 

Where {a PI ,a’ Pi }, {b Zi ,b ' zl } and {c Pi ,c* ; }are unknowns. 
Let 

(dx(0)p =Xpo^o(0 + yK,^(0 + >’p,^*(0} = Y T O(0 (18) 



(Bu(t)) p = i P Mt) + iU(t) + £M*(0 } = £ T d>(0 (19) 

(=1 

(£x(0)p = Wpo^(0 + ^{w„4(0 + w‘4* (0 } = 1C T <1)(0 (20) 

z=l 

Using equation (18)-(20) for each P,P = 1,2, 3,... p the 
right-hand side of equation (12) has the form 

ff T O(0 = (Y T +L t )O(0 

By equating the coefficients of the same- order Fourier 
series, we obtain for P = 1,2,3, . . .p and 

W T =(Y T +L J ) (21) 

Also, by using equation (1 1) we get the constraints 



F P0 — cipo x p (0) 



1 



1 ^1 






2^r J=l j 



L = 0 



( 22 ) 



And for i = 1,2,3, ... r we have 



F pi — a Pi + . c Pi — 0 

2 m 



in 



F Pi — a Pi | * _ c P0 + c Pi \L — 0 



(23) 

(24) 



3.2. The performance index approximation 

We now approximate the performance index J by using 



Fourier series. 

Let 

/7_ _ _ _* _* -*\T\ 



a- 



foo q, • ••A 4i 4 •••4) T 



a t’i ••• a ei a pi ■■■ Ar) j 

-k 4 4 -KP 



p= 



fe ^2 ^Zr bz . 1 ^2 • • • 



(25) 



(26) 



And define the matrices: 

®i(0 = 

<*> 2 (0 = 





\ 


H 




® T ( 0 y 

A 


V 


d> T (t) y 



(27) 



Note that a , /? and 0 , (£) and O 2 (^) are matrices of 
order ,P(2r + l)xl,Z(2r + l)xl and ,PxP(2r + l) and 
ZxZ(2r + 1) respectively. Using equation (15) and (16) 

and (25)-(27) the state and control vector can be 
expressed as 
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x(t) - 6 1 (t)a (28) 

«(0 = ® 2 (0£ (29) 



Substituting equation (28) and (29) in equation (14) we 
get 



7 = 2a T 6*(Z)S , d 1 (Z)a + 2a 1 



[f 



0^020,(0^ 



1 

+ — 

2 






® T 2 (t)R<S> 2 {t)dt 



(30) 



P 



Equation (30) can be computed more efficiently by using 
equation (10) and writing J as 

J=-a r [(^{L)^ c {L)®s\x+-a T {D®Q)a+-p T {D®R)P 

2! ^ 2 (3 1 ^ 

Where ® denotes the Kronecker product [25]. 



Remark: As it is shown in (31) the dynamic cost 
function is converted to a statistic cost function and it 
simplifies the problem considerably. 

3.3. Solution of the optimization problem 

The optimal control problem has been reduced to a 
parameter optimization problem which can be stated as 
follows. Find a and (3 so that J(a,/3) is minimized 
subject to the constraints in equations (21)-(24) 

Let 

L(a , p) = J(a, p) + £ X } F Pj (a, /?) + Z X) F*. (32) 

7=0 7=1 



'0 


f 


*1 




'1 


o' 


*1 


+ 


'O' 


_0 


0_ 


x 2 




_0 


1 


x 2 




1 



1 

x (o ) = _ yr 
2 . 

With the cost function 



7 = 



T I [x T (t)x(t) + u 2 (t) \lt 



(39) 

(40) 

(41) 



The problem is to find the optimal control u(t) which 

minimizes equation (41) subject to (39) and (40). We 
determine the direct Fourier series approximation for 
r = 10 
Let 

10 r 

(0 = «iofMO + ZK i<Pi(t) + } (42) 

/= 1 
10 , 

x 2 (t) = a 20 <j) 0 (t) + ZK4(0 + 44*(0 } (43) 

l=l 

10 , 

«(0 = bj 0 (t) + Z {brfii t) + Kf.it) } (44) 

/ = 1 

io , 

X x it) = c lo 0 o (t) + Z + Kfit) } (45) 

i = 1 
10 , 

i 2 (0 = c 20 ^o(0 + Z l c 2 /^( ? ) + c 2 ,y *(0 } (46) 

(=1 



Applying equation (30) we obtain the following 
approximation for J 



J — CL 1 rv + Cljr, + br 



4zk 

£ i = 1 



+ ^ 2 / + b?) + (afi 



-K)\ 



(47) 



where X-(\ ... \ X r ... /I*) represent the 

unknown Lagrange multiplier. Then necessary condition 
for stationary are given by 



dL 

da p 

dL 

da Pi 

dL 

db Zi 

dL 

w Zi 



^ (S' 

II 


+!>,• 

7=0 


dF Pj 

da p 


+ Z f 

7=1 


dF* Pj 

da p 


II 


7=1 


dF Pj 

da Pj 


7=1 


dF* p 

da Pj 


II 


+Z x j 

7=0 


dF Pj 

8b z , 


7=1 


dF* Pj 

~db~ Zi 


II 


7=0 


8F Pj 

db‘ Zi 


+ Z x j 

7=1 


dF* Pj 

db* z . 


= 0 


j = 


o,i,: 


>,... r 




= 0 


j = 


1,2 


,... r 





z = 0,1,2,... r P = 1,2 , . . .p (33) 

i = 0,1,2, . . r P=l,2,...p (34) 

z = 0,1,2 , ...r Z = 1,2 , . . ,z (35) 

i = 0,1,2, . . .r Z- 1,2 , ...z (3 6) 

(37) 

(38) 



Remark: As it shown in (8) the operational matrix P 
has so many zeros and this fact results in a considerably 
saving of computing time. 



4. Illustrative example 

Consider the linear singular system [26] 



Using equations (32)-(38) the direct method Fourier 
series approximation for x x {t),x 2 {t) and u(t) in 
equations (42)-(44) can be calculated. 

In table I and n , a comparison is made between the 
values of x t (^)and u(t) using present method 
with r = 10 , and the exact solution 



Table I .Estimated and exact values ofx x ( t ) . 



x \ (0 


Solution 


time 


exact 


Fourier series 


Number 




solution 


r=10 


1 


0.00 


1.00000 


0.94543 


2 


0.25 


0.70219 


0.69577 


3 


0.50 


0.49307 


0.49945 


4 


0.75 


0.34623 


0.34457 


5 


1.00 


0.24312 


0.24095 


6 


1.25 


0.17071 


0.17384 


7 


1.50 


0.11987 


0.11834 


8 


1.75 


0.08417 


0.08314 


9 


2.00 


0.05911 


0.06488 
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Table EL .Estimated and exact values of u(t) 





u 


(0 




Solution 


time 


exact 


Fourier series 


Number 




solution 


r=10 


1 


0.00 


0.70711 


0.66852 


2 


0.25 


0.49652 


0.49199 


3 


0.50 


0.34865 


0.35316 


4 


0.75 


0.24482 


0.24364 


5 


1.00 


0.17191 


0.17038 


6 


1.25 


0.12071 


0.12292 


7 


1.50 


0.08476 


0.08368 


8 


1.75 


0.05952 


0.05879 


9 


2.00 


0.04179 


0.04588 



5. Conclusion 

The Fourier series and the associated matrix of 
integration (operational matrix) are applied to solve the 
optimal control of singular systems. The method is based 
upon reducing a linear singular quadratic optimization 
problem to a set of algebraic equations. The advantage of 
using Fourier series as compared to Faguerre 
polynomials, shifted Fegendre polynomials and 
Chebyshev series is that in the present method the, most 
of the elements of the operational matrix given in (8) are 
zeros, and matrix D introduced in (10) is a diagonal 
matrix, hence making Fourier series computationally 
very attractive. As we can see from the tables, it is shown 
that with a relatively low r , we have a very accurate 
solution. 

6. Acknowledgements 

The authors wish to express their sincere thanks to the 
referee for valuable suggestions which improved 
the final manuscript. 

7. References 

[1] MS. Corrington. Solutions of differential and 
integral equations with Walsh function. IEEE 
Transactions on Circuit Theory, 20: 470-476, 1973. 

[2] C. Hwang and YP. Shih. Optimal control of delay 
systems via block pulse functions. Journal of 
Optimization Theory and Application, 45: 101-112, 
1985. 

[3] C. Hwang and YP. Shih. Solution of integral 
equations via Faguerre polynomials. Computers and 
Electrical Engineering , 9: 123-129, 1982. 

[4] M. Razzaghi and M. Razzaghi. Solution of linear 
two-point boundary value problem and optimal 
control of time-varying systems by shifted 
Chebyshev. Journal of the Franklin Institute, 
327:321-328, 1990. 

[5] HM. Pemg. Direct approach for the optimal control 
of linear time-delay systems via shifted Legendre 
polynomials. International Journal of Control, 43: 
1897-1904, 1986. 



[6] M. Razzaghi and M. Razzaghi. Fourier series direct 
method for variationl problem. Int. J. control, 48: 
887-895, 1988. 

[7] L. Dai, singular control system, Springer, Berlin, 
1989 

[8] M. Kuijper, First Order Representation of Linear 
systems, Birkhauser, Boston, 1994. 

[9] F. L. Lewis, A survey of linear singular systems, 
Circuits, Systems Signal Process. 5 (1): 3-36, 1986. 

[10] SL. Campbell. Singular Systems of Differential 
Equations. Pitman: London, 1980. 

[11] SL. Campbell, CD Meyer and NJ, Rose. Application 
of the Drazian inverse to linear systems of 
differential equations with singular constant 
coefficient. SIMA Journal of Applied Mathematics, 
31:425-611, 1976. 

[12] Z. Gao, and A. T. P, So. A general doubly coprime 
factorization for descriptor systems, Systems and Control 
Letters, 49(3): 213-224, 2003. 

[13] Gao, Z and D. W. C, Ho. Proportional multiple- 
integral observer design for descriptor systems with 
measurement output disturbances, IEE Proc. -Control 
Theory Appl, 151(3):.279-288, 2004. 

[14] Gao, Z. , and S. X, Ding,. Actuator fault robust estimation 
and fault-tolerant control for a class of nonlinear descriptor 
systems, Automatica, 43: (5), 912-920, 2007. 

[15] S.L. Campbell. On using orthogonal function with 
singular systems, proc. IEE Pt. D, 131: 267-268, 
1984. 

[16] F.L. Lewis and B. G. Mertzios , Analysis of singular 
systems using orthogonal functions, IEEE Trans. 
Auto Control, 32: 527-530, 1987. 

[17] W. Marszalek, On using orthogonal functions for 
the analysis of singular systems, proc. IEE Pt. D, 
132: 131-132, 1985. 

[18] P. N. Paraskevopoulos. Analysis of singular systems 
using orthogonal function. Proc. IEE Pt. D, 131: 37- 
38, 1984. 

[19] Z. Trzaska, Computation of the block pulse solution 
of singular systems, IEE Pt. D, 133: 191-192, 1986. 

[20] K. R. Palannisamy and Balachandran, Single-term 
Walsh series approach to singular systems, 
Int.J.Control, 46: 1931-1934, 1987. 

[21] K. R. Palannisamy and Balachandran, Analysis of 
time-varying singular systems via Single-term Walsh 
series approach proc. IEE Pt. D, 135: 461-462, 1988. 

[22] K. Balachandran, and K. Murugesan, Analysis of 
Different systems via Single-term Walsh series 
method, Int. J. Comp. Math. 33: 171-179, 1990. 

[23] K. Balachandran, and K. Murugesan, Analysis of 
electronic circuits using the Single-term Walsh series 
approach. Int. J. Electronics, 69: 327-332, 1990. 

[24] P. N. Paraskevopoulos, P. D. Sparis, and S.G 
Mouroutsos. Int. J. Systems Sci.16: 171-181, 1985. 

[25] P Lancaster. Theory of Matrices. Academic Press: 
New York, 1969. 

[26] K. Balachandran and K. Murugesan. optimal control 
of singular system via single-term Walsh series. 
International Journal of computer Mathematics. 43: 
153-159, 1992. 




23 




Automatic Control and System Engineering Journal, ISSN: 1687-4811, Volume 7, Issue 2, ICGST, Delaware, USA, Nov. 2007 



Ramazan Ebrahimi, He is a M.Sc 
student in control engineering in 
electrical Dept in Bahonar 
University of Kerman. 



Mahmoud Samavat, He got his 
B.Sc and M.Sc all in electrical 
engineering in 1982 and 1985 
respectively from North Dakota 
state university in USA. Also he 
got his Ph.D in control 
engineering from Dep of 
Automatic control and system 
engineering from university of 
Sheffield, England in Dec 2000. 
Since 1985 he is a faculty member 
of electrical Dept in Bahonar university of Kerman, his 
research area is stabilization, identification, model 
reduction and optimal control. 

Mohammed Ali Valim, He got 
his B.Sc in pure mathematics in 
1985 from Esfahan University in 
IRAN, he got his M.Sc in pure 
mathematics (Analysis) in 
Bahonar university of Kerman in 
1987. Also he got his Ph.D in 
pure mathematics (optimal 
control) in 1999 from Bahonar 
university of Kerman in IRAN. 
Since 1987 he has been a faculty member of Math Dep. 
Of Bahonar university of Kerman. His research area is 
analysis, Linear algebra and optimal control. 






24 





Automatic Control and System Engineering Journal, ISSN: 1687-4811, Volume 7, Issue 2, ICGST, Delaware, USA, Nov. 2007 




1308T 



www.icgst.com AOSI 



A Hybrid Genetic-Simulated Annealing Approach to TCPS based 
Hydrothermal System under Open Market System 



C. Srinivasa Rao 1 , S. Siva Nagaraju 2 , P. Sangameswara Raju 3 
1 Assistant Professor at G. Pulla Reddy Engineering College, Kumool, Andhra Pradesh, India 
2 Associate Professor at J.N.T.U College of Engg., Kakinada, Andhra Pradesh , India. 

3 Associate Professor at S.V.University, Tirupati, Andhra Pradesh, India. 

E-mail: csr_bit@yahoo.co.in 



Abstract 

This paper presents the analysis of Automatic generation 
control (AGC) of a two-area interconnected hydrothermal 
system with Thyristor controlled phase shifter (TCPS) in 
series with the tie-line under open market system. It also 
involves the optimization of integral controller 
employing a hybrid approach of modified Genetic 
algorithm and Simulated Annealing method (HMGASA) 
modified genetic algorithm method. Open transmission 
access and the evolving of more socialized companies for 
generation, transmission and distribution affects the 
formulation of AGC problem. So the traditional AGC 
two-area system is modified to take into account the 
effect of bilateral contracts on the dynamics. It is possible 
to stabilize the system frequency and tie-power 
oscillations by controlling the phase angle of TCPS 
which is expected to provide a new ancillary service for 
the future power systems. A control strategy using TCPS 
is proposed to provide active control of system frequency. 
Gain settings of the integral controllers without and with 
TCPS are optimized using the HMGASA technique 
following a step load disturbance in either of the areas. 
Analysis reveals that a TCPS is quite capable of 
suppressing the frequency and tie-power oscillations 
effectively as compared to that obtained without TCPS 

Key Words: AGC, Hydrothermal, TCPS, HMGASA algorithm, 
Disco participation matrix, open market system 

1. Introduction 

Large scale power systems are normally composed of 
control areas or regions representing coherent groups of 
generators. In a practically interconnected power system, 
the generation normally comprises of a mix of thermal, 
hydro, nuclear and gas power generation. However, 
owing to their high efficiency, nuclear plants are usually 
kept at base load close to their maximum output. Gas 
power generation is ideal for meeting the varying load 
demand. Thus the natural choice of AGC falls on either 
thermal or hydro units. A realistic situation may arise 
where an area regulated by hydro generation is 
interconnected to another area regulated by thermal 



generation. Nowadays worldwide, the electric power 
industry is in a transition from Vertically integrated 
utility (VIU) scenario, where a single utility owned and 
operated the generation, transmission and distribution 
systems and provided power at regulated rates to the 
deregulated scenario, where, competitive companies sell 
unbundled power at lower rates. Furthermore, various 
kinds of apparatus with large capacity and fast power 
consumption such as testing plants for nuclear fusion, 
steel factories etc increase significantly. When these 
loads are concentrated in power systems, they may cause 
a serious problem of frequency oscillations. Therefore, it 
is very important to consider how the control services of 
system frequency should be implemented. In a 
deregulated environment, any power system control such 
as frequency, will serve as an ancillary service. Thus, 
stabilization of frequency oscillations in an 
interconnected power system becomes challenging when 
implemented in the future competitive environment. A 
new frequency stabilization service which emphasis not 
only efficiency, reliability and economics, but also, 
advanced and improved controls for satisfying the 
requirements of power system operation, is much in 
demand. 

On the other hand, the concept of utilizing power 
electronic devices for power system control has been 
widely accepted in the form of Flexible AC Transmission 
Systems (FACTS) which provide more flexibility in 
power system operation and control [4]. This extra 
flexibility permits the independent adjustment of certain 
system variables such as power flows, which are not 
normally controllable. A Thyristor Controlled Phase 
Shifter (TCPS) is expected to be an effective apparatus 
for the tie-line power flow control of an interconnected 
power system In the analysis of an interconnected power 
system, some areas are considered as the channels of 
disturbances and in this situation, the conventional 
frequency control i.e., the governor may no longer be 
able to attenuate the large frequency oscillations due to 
its slow response [5, 6]. On the other hand, tie-line power 
flow control by a TCPS installed in series with the tie- 
line in between two areas of an interconnected power 
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system has the possibility to control the system frequency 
positively. The proposed control strategy will be a new 
ancillary service for the stabilization of frequency 
oscillations of an interconnected power system. 
Literature survey shows ample applications of TCP S for 
the improvement of dynamic and transient stabilities of 
power systems[3],[7]and[8]. However, no attempt has 
been made to improve the performance of AGC of 
interconnected power system under open market system 
considering TCPS in series with the tie-line. In the 
present work, analysis has been carried out for the AGC 
of interconnected hydrothermal power system under open 
market system considering a TCPS in series with the tie- 
line. An interconnected hydrothermal system involves 
widely different characteristics for the hydro and thermal 
subsystems. The characteristics of hydro turbines differ 
from steam turbines in that the relatively large water 
inertia used as a source of energy causes a considerably 
greater time lag in the response of the change in the 
prime mover torque to a change in gate position, and also 
an initial tendency for the torque to change in a direction 
opposite to that finally produced. The speed governor 
characteristics of the hydro units are widely different 
from that of the turbo governor. Moreover, the maximum 
permissible generation rate constraint for the hydro units 
is relatively much higher than the thermal units. Further, 
the effects of different Generation Rate Constraints 
(fairly slow response for the thermal plant and quite fast 
response for the hydro plant) on the selection of optimum 
controller settings for the thermal and hydro areas and on 
the system dynamic performance considering a TCPS in 
series with the tie-line is yet to be established. In view of 
the above, the main objectives of the present work are: 

1) To study the effect of a Thyristor Controlled 
Phase Shifter (TCPS) in a tie-line on the AGC 
dynamics of a hydrothermal system 

2) To develop a linearised model of a two 

area interconnected hydrothermal system under 
open market scenario considering a TCPS in series 
with the tie-line 

3) To optimize the gain settings of the integral 
controllers with and without considering TCPS 
employing HMGASA method 

4) To compare the dynamic responses with and 
without considering the TCPS. 

2.Tie Line Power flow model 
Considering TCPS 

Under open market system the power system structure 
changed in such a way that would allow the evolving of 
more specialized industries for generation (Genco), 
transmission (Transco) and distribution (Disco). A 
detailed study on the control of generation in deregulated 
power systems is given in [9]. The concept of 
independent system operator (ISO) as an unbiased 
coordinator to balance reliability with economics has also 
emerged [10,11]. The assessment of Automatic 
Generation control in an open market environment is 
given in detail in [12,13] and also provides a detailed 
review over this issue and explains how an AGC system 
could be simulated after deregulation. The recent 
advances in power electronics have led to the 



development of the Flexible Alternating Current 
Transmission Systems (FACTS). FACTS devices are 
designed to overcome the limitations of the present 
mechanically controlled power systems and enhance 
power system stability by using reliable and high-speed 
electronic devices. One of the promising FACTS devices 
is the Thyristor Controlled Phase Shifter (TCPS). A 
TCPS is a device that changes the relative phase angle 
between the system voltages. Therefore, the real power 
flow can be regulated to mitigate the frequency 
oscillations and enhance power system stability. In this 
study, a two-area hydrothermal power system 
interconnected by a tie-line under open market scenario is 
considered. 



Area 1 




2 Area 2 




h 



ZSi 



\n\z(si +</> 



\V2\^2 



Figure 1: Interconnected hydrothermal two area 

system with TCPS in series with the Tie-line 



Figure 1 shows the schematic of the two-area 
interconnected hydrothermal system considering a TCPS 
in series with the tie-line. TCPS is placed near area 1. 
Area 1 is the thermal area comprising of three reheat 
units and area 2 is the hydro area consisting of three 
hydro units. Without TCPS, the incremental tie-line 
power flow from area 1 to area 2 under open market 
system can be expressed as [1] 



2 ttT° 

^*12— r-w i-Att a) 

When a TCPS is placed in series with the tie line as in 
Fig 1, a current flowing from area 1 to area 2 can be 
written as 



hi 



|Ki|z(J 1 + ^)-|f 2 |z ^ 2 
jX 12 



( 2 ) 



From Figure 1 it can also be written as 



p tie 12 -jQtie 12 = \ v \\X-(Si + fi) 



|Fi|z(^ 1 + ^)-|f 2 |zj 2 
jX 12 



\ 



J 



Separating the real part of Eqn (3), we get 

\vAv 2 \ 

P tie\2 = ~ 7 “~L n (^l —S 2 +<!>) 

x \2 



(3) 

•(4) 



But in Eqn (4) moving S \ , 82 and ^ from their nominal 



values 8 ° , 8 ^ and (j)° respectively, we get 



^tiell 



nn 

Xu 



cos^ 0 -82 +^ 6> )sin(A S\ -A 82 +A (f>) 

(5) 

But ( A 8 \- A 82 + A(/>) is very small and hence, 



sin( A 8 \ - AS 2 + A(/>) « (A 8 \ - A 82 + A <p) 
So Eqn (5) can be written as 
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^tie 12 — 



NN 

*12 



cos(^° -8° + <j>° )(A(J| -AA 2 + A^) 



( 6 ) 

Let L| 2 = cos(A 2 ° - + / ) (7) 

Thus,Eqn.(6) reduces to kPtiell = ^12 (^^1 + A#>) 

( 8 ) 

Therefore, AP^^ = ^12 (AAi - AA2 ) + T^A(p (9) 

We also know that, A8\-2n\Af\dt and 

A8 1 =2n\Af 1 dt (10) 

From Eqns.(9) and (10), it can be written as 

AP tfc i2 = 2^2 (j Af x dt - J Af 2 dt) + P 12 A^ (11) 

Taking the laplace transform of Eqn.(ll), it can be 
written as 



A/W12 (.s) = [AFj (5)- A F 2 (s)]+ T 12 A<p(s) (12) 

As per Eqn.(12), it can be observed that the tie-line 
power flow can be controlled by controlling the phase 
shifter angle A cp . The phase shifter angle A cp (s) can be 
represented as [12-14]: 



A<p(s) ■■ 



K 



* 



l+^r 



- A Error\ (s) 



•(13) 



ps 



Where K ^ and T ps are the gain and time constants of the 
TCPS. Thus, Eqn.(12)can be rewritten as 

AP lie l 2(s) = — [A/) (s) - AF 2 (s)] + 7) 2 AErro^s) 



1 + sT, 



ps 



(14) 

If the frequency deviation A f x is sensed, it can be used 
as the control signal (i.e., A Error x - A f \ ) to the TCPS 

unit to control the TCPS phase shifter angle which in turn, 
controls the tie-line power flow. Thus, 

Aq*s) = - AF!(j) (15) 

l + sT ps 

and the tie-line power flow perturbation becomes 



APfieU (J) = [Afj (5) - A F 2 (s)\+T [2 A F x ( 5 ) 

S l + sT ps 

(16) 

But in the open market system the actual tie-line power 
flow also includes the demand from Discos in one area to 
Gencos in another area. It can be represented as follows. 
A P t ie 1-2, actual = &Pfie \2 ( s ) + (demand of Discos in area 2 
from Gencos in area 1) - (demand of Discos in area 1 
from Gencos in area2) 



3. System Investigation 

The AGC system investigated is composed of an 
interconnection of two areas under open market system. 
Area 1 comprises of a reheat system and area 2 comprises 
of hydro system. The detailed transfer function models of 
speed governors and turbines are discussed and 
developed in the IEEE committee report on Dynamic 
models for Steam and Hydro turbines in Power 



systems [14]. The detailed small perturbation transfer 
function block diagram model of two area hydrothermal 
system under open market scenario along with the 
incremental model of TCPS in series with the tie-line is 
shown in Figure 2. Nominal parameters of the TCPS 
system are given in the Appendix. 

4. Optimization of Integral Gain Settings 

Modified genetic algorithm method is used to obtain the 
optimum integral gain settings. A performance index 

which is denoted by J = J + AP^^dt 

0 

is minimized in the presence of GRCs to obtain the 
optimum values of Kf\ and Kf 2 . A value of 0.65 is 
considered for both a and /? . The modified genetic 

algorithm method resembles the normal traditional 
genetic algorithm which takes into account the nature of 
genetics [15]. But in modified genetic algorithm there are 
3 modifications i.e (a) Modification in parent selection (b) 
Improvement in child string (c) Elitism. The algorithm 
for modified genetic algorithm method can be written as 
shown below. 

Step l:As the genetic algorithm takes pairs of strings, we 
create a random number of strings depending 
upon our necessity and also note down their 
decoded values along with setting a maximum 
allowable generation number t max 
Step 2: Using the mapping rule we next find out the 
corresponding values of Kj for 
the corresponding decoded values. 

Step 3:Using these values of Kj the fitness function 
values are found out. 

Step 4: Next the process of reproduction is carried out 
on the strings to create a mating pool 
Step 5: The process of crossover and mutation is 
also carried out on the strings with probabilities 
of 0.8 and 0.05 respectively. 

Step 6: After the termination criteria is met with, 
the value of Kj with minimum fitness 

function value is considered as optimum value. 
Here the following modifications are applied for the 
selection of parents so that the strings with large values 
of fitness are copied more into the mating pool. 

1 . The first parent in each reproduction is the string 
having the best fitness value. The second parent is 
selected from the ordered population using 
normal selection technique. 

2. At the i th reproduction, first parent is the i th string of 
the population arranged in the order of fitness 
values. Second parent is selected from the ordered 
population using normal selection technique. 

To improve the quality of child strings, the strings having 
fitness values less than average fitness values of the 
previous generation are identified first. Their quality is 
improved by selecting a number of bits at random with 
probability P m i and complimenting. In order to avoid 
losing of best results during the selection, elitism is used. 
So here the best solution of every generation is copied to 
the next so that the possibility of its destruction is 
eliminated. 
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Figure 2: Two area hydrothermal system under open market scenario 



A new hybrid method based on the combination of 
modified Genetic algorithm and simulated annealing 
has been proposed wherein the simulated annealing 
takes over the optimization of integral controller once 
the genetic algorithm terminates. This new method is 
called as hybrid modified Genetic- Simulated 
annealing method (HMGASA) which indeed gives 
good results than the techniques considered separately. 
The simulated annealing method resembles the cooling 
process of molten metals through annealing. The 
algorithm for simulated annealing method can be 
written as shown below. 



Step 1: Choose an initial point x^ , a termination 
criterion s .Set T a sufficiently high value, 
number of iterations to be performed at 
a particular temperature n ,and set t=0. 

Step 2: Calculate a neighbouring point 

x 0+l) =7V(x^) . Usually a random point in 
the neighborhood is created 
Step 3: If A£ = £(^ +1 ) )-£(j W )< 0, set t = t+1; 



else create a random number (r) in the range 
AE 

(0,1). If r < exp(- ~Y~)set = t + 1; else go to 



step 2 



Step 4: If 




< £ 



and T is small, 



Terminate else if (t mod n)= 0 then lower 
T according to a cooling schedule. And go to 
step 2. Else go to step 2. 

Optimum values of integral gain settings of area 1 and 
area 2 without and with TCP S in series with the tie-line 
are tabulated in Table- 1 



T abl e 1. Optimum values of Integral Gain Settings 



Area 


Without 

TCPS 


With TCPS 
inseries 
with tie-line 


Thermal 


kfi = 0.650 


k a =0.940 


Hydro 


£,2=0.470 


£,•2=0.570 



From Table-1 it is seen that the optimum values of the 
integral gain settings in areas 1 and 2 considering 
TCPS are higher than those obtained without TCPS. 



5. Dynamic Responses and Discussions 

Considering GRCs, simulation studies are performed to 
investigate the performance of the two-area 
hydrothermal system under open market system 
without and with TCPS in series with the tie-line. A 
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step load disturbance of 0.04 pu MW is considered in 
either of the areas. Also an additional case is also 
considered when contract violation occurs in either 
area. So in this contact violation an additional load of 
0.03 pu MW is considered in both the areas after the 
time span of 30 sec. In this contract violation the 
Gencos which are present in that particular area where 
the violation has taken place indeed take up that extra 
load while the remaining generators of other area 
generate the power which they had been generating 
before. It can also be seen that this contract violation 
does not effect the systems frequency response and 
altogether the system comes back to its original state 
after some disturbance initially. 




time (s) 

Figure 3: Variations in area frequencies and tie-line 



power 

It is seen from Figure 3 that, with the TCPS in series 
with the tie-line, dynamic responses for A/q , A /2 have 

improved significantly and it can be observed that the 
oscillations in area frequencies have decreased to a 
considerable extent with the use of TCPS. Depicted in 
Figure 4 are the generations of Gencos in first area and 
Figure 5 depicts the generation of Gencos in second 
area for a load of 0.004 p.u Mw in both areas. This idea 
has been implemented in the block diagram and it has 
also been represented in block diagram of Figure 2. 



x 10' 3 





x 10' 3 




Time (s) 

Figure 5: Generations of Gencos in area 2 



Figures 6-8 shown below show the case of contract 
violation where an additional load of 0.03 pu MW is 
demanded in each area. As mentioned earlier it can be 
seen in Figure 6 that the frequency error has some 
transients and after some time it comes back to its 
original zero position after some time. Figures 6-7 
depict the generations of Gencos in both areas during 
contract violation and also with the use of TCPS the 
oscillations are reduced drastically 





Figure 6: Variations in area frequencies and tie-line 
power deviations during contract violation 




Figure 7: Generations of Gencos in area 1 during 
contract violation in area 1 
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Figure 8: Generation of Gencos in area 2 during 
contract violation in area 2 



Depicted in Fig 9 are the frequency deviations for the 
various values of integral controllers generated during 
the process of HMGASA algorithm. The optimum 
values of integral controllers without and with TCPS 
have already been tabulated in Table 1. Here also the 
oscillations in area frequencies have been greatly 
reduced by the use of TCPS 




various values of integral controller 



6. Conclusions 

In this paper, a tie-line power flow control technique 
by TCPS has been proposed for a two-area 
interconnected hydrothermal power system under open 
market scenario. Gain settings of the integral 
controllers are optimized using HMGASA algorithm 
in the presence of GRCs by minimizing a quadratic 
performance index. A control strategy has been 
proposed to control the TCPS phase angle which in 
turn controls the inter-area tie-line power flow. 
Simulation results reveal that frequencies and tie- 
power oscillations following sudden load disturbance 
in either of the areas can be suppressed by controlling 
the phase angle of TCPS. It may be therefore 
concluded that, the tie-line power flow control by a 
TCPS can be expected to be utilized as a new ancillary 
service for stabilization of frequencies and tie power 



oscillations in the deregulated environment of power 

systems. 

7. Appendix 

All the notations carry the usual meanings 

TCPS data: 

Tp S 

K ^ - 1.5 rad / Hz 

^max ~ 

Anin - _ 10 
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Abstract 

This paper deals with the design of a dynamic output 
feedback controller for nonlinear and uncertain 
systems. The proposed controller is based on the 
Takagi-Sugeno (TS) fuzzy approach. Both stabilities 
in a pre-specified LMI stability region and Hqo 
control synthesis have been developed. The control 
laws are derived considering the so-called dynamic 
output feedback parallel-distributed compensation 
(DOPDC). The resolution of the problem associated 
with the controller synthesis has been achieved 
considering an LMI approach. The capability of the 
proposed approach is illustrated through a numerical 
example. 

Keywords: Fuzzy system , Robust control, dy- 
namic output feedback, robust D-stability, linear ma- 
trix inequality (LMI). 



1 Introduction 

Over the past two decades, fuzzy logic control has 
been given an increasing attention. The implementa- 
tion of fuzzy control strategies has led to successful 
applications. Of particular interest are the Takagi- 
Sugeno (TS) fuzzy models which are currently con- 
sidered as universal approximator for any nonlinear 
function [22]. The (TS) fuzzy model represents a non- 
linear system by a set of local linear models, which is 
smoothly blended together through membership func- 
tions. This representation offers the possibility of us- 
ing conventional linear systems for the analysis and 
design of the fuzzy controller. Being an interesting 
procedure of the control design framework, the par- 
allel distributed compensation (PDC) is applied in 
the present paper. The PDC control structure uses 
a nonlinear state feedback controller, which mirrors 
the structure of the associated (TS) model. The gains 
of the controller can be determined automatically us- 



ing an LMI formulation [1, 3, 5, 6, 13, 20, 23]. Since 
the parametric uncertainty is a crucial factor from 
the point of view of stability and performance, it’s 
very important to consider the robust stability against 
parametric uncertainty in the (TS) fuzzy controller 
[4, 7, 10, 11, 21, 23, 26]. Recently, many works on 
(TS) fuzzy model-based control of nonlinear systems 
were developed to stabilize the (TS) fuzzy model with 
an adequate level of performance such as the H 2 per- 
formance [8], Hoc performance [10, 13, 14, 17, 20] or 
the mixed H^/Hqo performance [9]. In this paper, we 
extend the LMI based approach to develop a dynamic 
output feedback controller for a class of fuzzy uncer- 
tain systems. These are described by (TS) fuzzy mod- 
els in which both robust D-stability and prescribed 
Hoo performance are required to be achieved. The 
paper is organized as follow: 

Problem formulation and definitions are presented in 
section 2. The analysis and the synthesis of the re- 
sults of the proposed fuzzy dynamic output feedback 
controller, which guarantees the 11^ performance and 
quadratic D-stability, are developed in section 3. In 
section 4, numerical example is provided to confirm 
the design effectiveness and to illustrate the desired 
performance. 



2 System description and defini- 
tions 



A convenient and flexible tool for handling complex 
nonlinear systems is the (TS) fuzzy model, where the 
consequent parts are linear models connected by IF- 
THEN rules [3, 5, 6]. A discrete-time (TS) fuzzy 
model with uncertainties is composed of r plant rules 
that can be represented as 
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Plant Rule i: 



The overall (DOPDC) controller can be rewritten as 



If 0i(k) is Mi i and . . . and 0 p (k) is M Pi Then (1) 
x(k + 1) = (Ai + AAi)x(k) + (Bi i + AB\f)u(k) 

+ (_B 2i + AB2 i )w{U) 
z(k ) = + Di.u(k) 

y(k) = C 2i x(k) + D 2i w(k) 
for i = 1, 2 , . . . , r 

where Mj i is the fuzzy set; x(k) G IR n is the state, 
u(k) G IR m is the control input, y(k) G JR P is 
the output, z(k) G JR q is the control output and 
w(k) G IR* is the exogenous disturbance input with 
{w(k)} G L 2 [ 0, oo). 

A i} B i., B 2i , Ci i ^Di i ,C 2 i and £> 2 . are known real con- 
stant matrices and A A^ AB\. and AB 2 i are time- 
varying but norm bounded uncertainties with the fol- 
lowing form: 



x(k + 1) = £ hj(0(k)) Aj x(k) + Bj y(k) ( 6 ) 

3 = 1 ' ' 

r 

u(k) = Y J h j m))c j m 

3 = 1 

The closed-loop system obtained by connecting the 
dynamic controller ( 6 ) and the uncertain fuzzy system 
(3) can be expressed as 

r r 

x(k + l) = hih i 

i=i j = 1 
r r 

hihjCij x{k) 

i=i j = l 

where 



A ii7 - x(fe) + Bij w(k) (7) 



[AA< AB l4 AB 2i ] = HF{k) [E t E u E 2i ] ( 2 ) -j 



where H, Ei , Bd. andB 1 ^ are known constant matrices 
of appropriate dimensions and F(fc) is an unknown 
real time varying matrix with Lebesgue-measurable 
elements satisfying F T (k)F(k) < I. 

The dynamic fuzzy model can be represented by the 
following overall model, which combines all the local 
models through membership functions. 



*<* + !) = £ h i {e(k))UA. l + AA i )x{k) 

i = 1 ' 

-\-(Bi i ~\-ABi i )u(k) + (.B 2i + AB 2 i )w[k) 

z (k) = h i( 0 ( k )) ( c U x ( k ) + Di^fc) ] ( 3 ) 

2=1 ' ' 

2/(fe) = X] h i( 0 ( k )) ( C 2 iX(k) + D 2i w(k) 



2=1 



A-ij — + HF(k)Eij , Bbj — + HF{k)E \ i , Cdj — CV 7 

(8) 



with 



^- 2 J — 



B jC 2 i Aj 



, -Bij — 



BjD 2i 



Cij = 



5 


II 

ie? 


EiE u C, 


5 


0 


, and 


& I 

F? 

II 


H) H 



C u D u C 3 
E u =E 2i , = 



Now, we recall the following definitions. 

Definition 2.1 : Suppose ^ is a given positive num- 
ber. A system of the form (7) is said to have B 2 gain 
less than equal to 7 if 



k = 00 



k=oo 



^ z T (k)z(k) < y 2 w T (k)w(k) 



(9) 



fc =0 



fc =0 



where hi(6(k)) is the normalized weight for each rule, 
that is: 

r 

0 < hi(6{k)) < 1 and hj{6{k)) = 1 (4) 

2=1 

for i = 1 , 2 , • • • , r 

Conformably to the description ( 1 ), it is quite nat- 
ural to seek a dynamic output feedback parallel- 
distributed compensation (DOPDC), [ 6 , 24, 25] in the 
form: 

If 0i(k) is M\. and . . . and 0 p (k) is M p . Then (5) 
x(k + 1 ) = Aj x{k) + Bj y{k ) 
u{k) = Cj x{k) 
forj = 1 , 2 , ...,r 



Definition 2.2 ; A subset D of the complex plane 
is called an LMI region if there exists a symmetric 
matrix Y G 1R 9X9 and a matrix II G 1R 9X9 such that 

V = {z — x + iy G (7; fx>(z) < 0} (10) 

where the characteristic function fx> is given as fol- 
lows: 

Mz) = r + nz + n T z (11) 

The uncertain system x(k-\- 1 ) = (A-\- AA)x(k) is said 
to be quadratically D-stable if there exists a symmetric 
positive- definite matrix P such that for all admissible 
perturbation AA , the following matrix inequality 

r <g> P + n <g> P(A + A A) + U T 0 (A + A Afp < 0 

is true, where the LMI region is defined in (11) 
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Problem formulation: The aim of this paper is to de- 
sign a full order fuzzy dynamic output feedback con- 
troller of the form ( 6 ) such that the resulting closed 
loop fuzzy system (7) is quadratically D-stable in the 
given LMI region while satisfying a prescribed 
performance level irrespective of the parameters un- 
certainties. 



(I*® Mas) -e 0ij I* | <0 (15) 

\(Li 0 HjP) 0 0 -£7)1 j 



1 < i < j < r 
where 



3 Dynamic output feedback 
fuzzy controller design with 
D-stability constraints 

In this section, we develop the (DOPDC) controller 
( 6 ) for uncertain (TS) fuzzy model (3), taking into 
account the formulation of the previous section. 

We first give the following results which will be used 
in the proof of our results. 

Lemma 3.1 [25] Given any matrices X , Y and Z 
with appropriate dimensions such that Y > 0 then, 
we have 

X T Z + Z T X < x t yx + z t y~ x z 

Lemma 3.2 [15, 26] Let A , H , E , P and F be real 
matrices of appropriate dimensions such that P > 0 
and F t F < I . Then, for any scalar £ > 0 such that 
P — sHH t > 0, we have 



©ii = r (8) p + n (g) PAa + (n <g> pa u ) t , 

(pij = 2 r ® p + n ® P{Aij + Aji) + (n ® p(Aij + Aji ))^ 



Hc = 



H H 



5 M c ij — 



Eij 

E 






and II = Lj L 2 



Proof: Firstly, we show that for any non zero w(k) E 
L 2 , the uncertain system (7) satisfies (9) under zero 
initial condition. To this end, we introduce the crite- 
rion: 



J = ^ ( z T (k)z(k ) — 7 2 w T (k)w(k)) (16) 

k = 0 



Noting the zero initial condition and the system (7), 
we can deduce 



(A + HFE) t P-\A + HFE) <A T (P - eHH T )~ l A J = £ [ V (x(k + 1)) - V(x(k)) + z T (k)z(k ) 



+ e~ 1 E T E 

Theorem 3.1 For a prescribed performance 7 > 
0 and LMI stability region V, the system (7) is D- 
stable with H ^ norm bound 7 if there exists a matrix 
P > 0 such that satisfying the following matrix in- 



fc =0 



— 7 2 w T (k)w(k)] —V(x(oo)) 



< 



^ [A V(x(k)) + z T (k)z{k) — 7 2 w T (k)w(k)] 



k = 0 

00 s r r r r 



equalities: 










/ — P * * * 


* 


* ^ 






0 — 7 2 I * * 


* 


* 






An Bu -P _1 * 


* 


* 


< 0 




Cu 00 -I 


* 


* 




IP 

0 

0 


—^iil 


* 






^ 0 0 H T 0 


0 


-epl) 






/ e« 


* 


* \ 






(L2<S>Eii) —Soul 


* 


< 0 




\(ii 0 # T P) 


0 - 


-n l i) 






i = l,...,r 








- 


-4 P * * 


* 


* 


* 




0 — 47 2 / * 


* 


* 


* 


Aij 


+ Pzj + Pjz — P 


1 * 


* 


* 


Cij 


T 0 0 


-I 


* 


* 


Eij E u 0 


0 


—£ijl 


* 


Eji Ei 0 


0 


0 - 


£ijl 




0 0 Hj 


0 


0 


0 



— s y^ hihjh u h v (fAijx{k) + 

/c=0 ^ i=l j=l ix=l i>=l 

x P[^x(fc) + B uv w(k)\ + [ Cijx(k)} T [C uv x(k )]) 



( 12 ) 



x 1 (k)Px(k) — 7 2 w 1 ( k)w(k ) 



(17) 



where V(x(k)) = x T (k)Px(k) is a Lyapunov function 
candidate. 

By Lemma (3.1), we have 



* 

* 

* 

* 

* 

* 

~-i 1 



(13) 

\ 



r r r r 



< 0 



hih,h u h v [ C ij x(k)\ T [C uv x(k )] 

i=l j = l n=l p=l 
r r r r 

=4 E E E E hihjKKiw^ + 

i= 1 j = l n=l p=l 

x[(C IIP + C„„)x(fc)] 

1 r r 

<TE hihj[(Cij + C ji )x(fc)] T [(C ij + Cji)x(fc)] 



(14) 



»=i 1=1 



(18) 
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and 



r r r r 



Then, applying the Schur complement to (14) gives 
P- 1 - SijH c Hj > 0 (24) 

and 



(r ciJ -+r ci ,) (p SijH c H c ) (r c ^+r c j^)-b£- E ci jE clJ 



EEEE hihjh u h v [Aijx(k) Bijw{k^ 

i= 1 j el n=l p=l 

x P[ 7 f ni; x(/c) T P wv ie(^)] 

^ r r r r 

~ ^ ^ ^ ^ ^ ^ ^ ^ hihjh u h v [(Aij T Aj^x(E) 

i= 1 j=l n=l p=l 

+ (^ij H - Pji)ie(/c)] T 

x P[{A UV + A vu )x(k ) + ( B uv + P vu )ie(/c)] 

^£X>; [(t 4 ^ + Aji)x(k) + (Pij + Pji)re(/c)] T (Y C 2jpr c ji+P c P(A;)P c ij) P(r c ^+r c j^+P c P(A:)P c ^j) 
2 = 1 j = l 

x P[(A^ + Aji)x(k) + (Bij + Pj?i)ie(/c)] 



+ (Y c ij + T cj ^) t (T c ^ + Y C ji) + 4P 7 < 0 (25) 

Considering (24) and using lemma (3.2), we have for 
1 < i < j < r 



(19) 



+ (Y C ij + Y c ^) T (Y C ij + T CJ -i) + 4P 7 
= (r cij + ^cji) T P(Vcij + T cji ) + 

(Y cij + Y cji) T (Y cij + Y cji) + 4P 7 < 0 (26) 



Therefore 



oo f r r 

J - E 4 { E E hih i ( 



where F(k) 



A — ( ) V 2=1 J = 1 

[{Aij + Aji)x(k) + {Bij + Pji)re(A;)] T 
x Ppy + + (Bij + Bji)w(k )] ( 2 °) 

[(C^j + Cji)x(/c)] T [(Cij + Cji)x(/c)] 



'P(Jfe) 0 

0 F(k)_ 

Using (12), and following a similar way to the deriva- 
tion of (26), we can obtain 



r T pr + t t Y + p <r o 

1 C22 r i C22 T -L -L C 22 W -*7 ^ U 



(27) 



— 4 x T (k)Px(k) — 4y 2 w T (k)w(k) 
This implies 



Therefore, the inequality in (21) together with (26) 
and (27) gives that for any £(fc) ^ 0 J < 0, which 
implies X)fc=o° zT (k) z (k) < l 2 Y%Z^ wT {F)w{k) f° r 
any zeros w{k) E P2[0,oo) Next, we establish the D- 
stability performance, we assume w{k) = 0. Then the 

J < j E ( t ± hAfm «%,J + T eji ) I '(T« + T eJ( f 0Sed - l00 > ) 8y8tem (7) is S ‘™ by 

r r 

x(k + 1) = hihjAij x{k) (28) 

2=1 J = 1 

Now, we introduce 



k = 0 v 2=1 j=l 
\T j 



+ (r C 2j + r C ji) T p(r C ij + r C ji) + 4P 7 )£(fc) 

oo ( r 

= E E h l C^cK + r^iTeii + P 

fc=0 ^ i=l 



+ 2 



oo ^ r r 

^e{ee 

k = 0 ^ 2=1 j>i 



Vx i x \T/^ i x \ ^mj — hi hj ( T(8)P+II(8)P^4ij+(n(8)P7lij)" 



hihjt ■ (fc) (Y c ij + Y c ji) (Y c ij + Y c ji) 



+ (r C ij + r C ji) T p(r C ij + r CJ -i) + 4 p 7 ^(fc) 

where 

Y C 2j = [C^ 0] , T c ij = [ A^ Pij] , P 7 = — 
and £ T (^) = [x T (&) (fc)] 

Set for 1 < i < j < r 



} 



2 = 1 j = l 

r r 



2 = 1 j = l 



— — ^ ^ ^ ^ hihj i 2 T (8) P + n 0 P{Aij + Aji) 



( 21 ) 



P 0 
0 7 2 / 



+ (n ® p(Aij + Ajj)) r 



= ^/i 2 (r0P + n0 p^i + (n 0 pa ^)) 5 



r r 



( 22 ) EE hihj l 2 T(g)P +II®P(Aij T Aji)+(Il 0 P(Aij T ^ 4 ji))^ 

2=1 j>i ^ 



(29) 



T — 

-L CIJ — 



CijO 



r — 

5 1 C2J 



Tlij Bij 



•> Ndj — 



Eij Fiij 



and E c ij = 



cij 



CJI 



On the other hand by the Schur complement, we have 
from (15) 



(23) 



(t> i j+e 0 (jL2®M cij ) T {L 2 ®M cij )+eo i jL(®PH c )(Ll®PH c ) T < 0 

(30) 
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Using the matrix fact HFE + E T F T H T < eHH t + where 
£~ 1 E T E, we have 



Aj + (Lj 0 PH C )(1 0 F{k)){L 2 0 M cij ) 

+ (L 2 0 M c ij) T (I 0 F(k)) T {Lj 0 PP C ) T 
= c hj+Li L 2 ®PH c F(k)M cij +(Lj L 2 0PH c F(k)M cij ) T 
= 2r 0 P + n 0 P(Aij + -f H c F(k)M c ij) 

+ (II (8) P(Aij + + H c F(k)M c ij)) T 

= 2T0P +II(8)P(Aij + Aji)+(Il(8)P(^4.ij + Aji)) T < 0 

(31) 

Using (13), and following a similar way to the deriva- 
tion of (31), we can obtain 

r ^ p + n ^ + (n 0 pa^)) t < o (32) 

Then, it follows from (31), (32) and (29) that 
F^.. < 0. Hence, the uncertain system (28) is 
quadratically D-stable. This end the proof. □ 

Now, we are in position to give the main re- 
sults on the solvability of the robust output feedback 
Hoc control problem with pole placement constraints. 
Theorem 3.2 Consider the uncertain discrete fuzzy 
system (3). Given a scalar 

7 > 0 and any LMI stability region V, then there ex- 
ists a full- order fuzzy dynamic output feedback con- 
troller (6) such that the closed-loop system (7) is 
quadratically D-stable and the prescribed H 0 c perfor- 
mance is guaranteed if there exist matrices X > 0, 

Y > 0, and \Eq, 1 < i < r, such the fol- 

lowing LMIs hold for some scalars Eij and £q • 



( J i ^ ^ ^ ^ ^ 

0 — J2 * * * * 

Mi.. M hii —J\ * * * 

M2 0 0 —I * * 

^zz 

m 3u M 6i 0 0 -J 3ii * 

0 0 N 0 0 -JiJ 



< 0 



(33) 



$ii * 

L 2 0 M 3ii — Jr >i 
Li 0 N 0 




<0 



(34) 



/— 4 J 


l * * 


* 


* 


* 


* 


* 


* \ 


0 


— 4J2 * 


>!< 


* 


* 


* 


* 


* 


1—1 1 ij 


S 5y -X 


>!< 


* 


* 


* 


* 


* 




0 0 


— I 


* 


* 


* 


* 


* 


M^ 


; M 6i 0 


0 - 


J 3 ij 


* 


* 


* 


* 


Mz it 


; Me ] 0 


0 


0 - 




* 


* 


* 


0 


0 M 


0 


0 


0 


-Ji 


.. . * 


* 


M 4iJ 


O 

O 


0 


0 


0 


0 


J -I 


* 


0 


0 Mh. 

I ZJ 


0 


0 


0 


0 


0 


-d 




7 <Pij 


* 


* 


* 


*\ 








L 2 0 Mgy 


- J 5io 


* 


* 


* 








Li 0 M 


0 




* 


* 


< 0 






L 2 0 M 4 y 


0 


0 


-I 


* 








\Li ® Ml. 


0 


0 


0 


-d 







< 0 



(35) 

(36) 



Ji = 
J. 



Sly - M Uj + Mly, S 2 y - M 2 y P M^, ^y ~ y A M^, 

<Pij — 2 r 0 j\ + n 0 (Afiy + iiii^) + n T 0 

$u — r 0 </i p n 0 ii^i ii + iV 0 M 4ii 

/] „ , 

X ’ ^2 1 h Fsij £ijh i 6^ I 

r 4 y. =diag(£j>I,£j>I), J 5i =e 0ii I, J 5ij = diag(s 0ij I,£ 0ij I), 
■hi = £oy h Jeij = diag{e^I,£^.I) 

_ AY + Bi'S'j Ai 

ly [ A XA, + t> :j C 2l 

M 3 y = [EiY + ^VjEi], M 4 y = 



M 2 y = [Ci + A^-Ci] , 



P, - 0 

(C 2i -C 2 .)F0 



M 5 .. 

Ozj 

m 7 .. 

1 ZJ 

N = [H T H t X } , M = 



D>. 

XB 2i + i’jDi, 

0 0 

X(B U - B h ) <l> :l - <h 



, M 6i = [E 2i ] , 

, M 8i . = 



' Ms , 
M 3 ,' 



'H t H T X 

h t h t x 



Furthermore, a desired robust dynamic output feed- 
back Hqq controller is given in the form of (6) with 
parameters as follows: 



A = S~\A - XAiY - XB u ^i - <S>iC 2i Y)W- T 

(37) 

Bi = S- 1 ^ (38) 

Ci = W T (39) 

where S and W are any nonsingular matrices satisfy- 
ing 



SW 1 = I — XY 



(40) 



Proof: Recall that our goal is to derive the expressions 
of the controller parameters from (12)- (15). To do 
this, we partition P and P _1 as 



<"> 



where the partitioning of P and P -1 is compatible 
with that of Aij defined in (8). 

Now define 






( 42 ) 
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which imply that PT\ = T 2 and T ^ PT\ = T f T 2 
By applying the Schur complement to (35), we obtain 



( — 4«/i * * * 

0 — 4J2 * * 

S ly S 5y ~Jl * 

s 2i . 0 0 — / 

M Zij M 6ij 0 0 - 

0 0 M 0 

T 



0 \ 

0 

m 7 .. 

1 zj 

0 

0 

V 0 / 



/ 0 \ 

0 

m 7 . 

1 1 

0 

0 

0 ) 



* \ 

* 

* 

* 

-J3.. * 

O t j 

0 “ J 4 ij 

( M l\ 

0 
0 
0 
0 

V 0 / 



Now, by applying the Schur complement to (36), we 
obtain 

<Pij * 

L 2 0 M Sij -Jbij 

Li0 M 0 - 



(47) 



( M L\ 
0 
0 
0 
0 

V o / 



< 0 




(43) 



By lemma (3.1), we can be deduce from (43) that for 
1 < i < j < r 



<0 



Therefore, by lemma (3.1) and by performing a Kro- 
necker product property, the inequality (47) implies 
that 



<Pij * 

1/2 (8) — J5 

Li <g> M 0 



( — 4Ji * * * 

0 — 4J2 * * 

a + M-ij H 5y - Ji * 



'2 zj 

Ms, 

0 



0 

Me, 

0 



* 

* 

* 

0 -I * 

0 0 -J 3 . 



\ 



<0 (44) 



M 0 0 -J Ai J 



where 



A4ij — 




(48) 



* 

1/2 0 Mg.. J$ij * 

L 1 0M 0 -J 6 ... 



0 



0 



X(B U - - T,) + (^ - ^)(C 2i - C 2j )Y 0 

= M 7 . . m 4 . . 



where 

2T ( 



1 Ji +11(8) (^iij + Mtij) + (II (8) (^1 ^ + Mij)) 1 



We can verify that the matrices inequalities in (48) 
and (34) respectively can be written as 

diag((I (8) T^), /, I) x (15) x diag((I (8) if), /, I) < 0 

(49) 

diag{{I (8) if ), /, /) x (13) x diag((I (g) if), /, /) < 0 

(50) 

From (45), (46), (49) and (50), we conclude that 
with controller parameters (37)- (39) the closed-loop 
system (7) is robustly quadratically D-stable and (9) 
is satisfied. This end the proof. □ 

Remarks 3.1 1. Theorem (3.2) provides a suffi- 

cient conditions for the solvability of the robust 
Hqq output feedback control problem with pole 

diag(T f, J, T 2 t , I, I, I, I) x (14) x diag^I, T 2 , J, J, 1,1) < 0 P^ement constraints. We note that (33)-(36) 
1 2 y V ’ yV ' are LMIs in X > 0, Y > 0, and 

\Eq, 1 < i < r, when 7 > 0, > 0 

are given. In this case, these LMIs can be solved 
efficiently by resorting to some standard numer- 
ical algorithm. 



Using the relationship 

(XAiY + XBuVi +*&<¥) 

+ ( XAjY + XB^t + i>,C 2 ,Y) = 

(XAiY + XB li * i + § i C 2i Y) 

+ ( XA,Y + XB Li * :) + <f>jC 2j Y)+ 

X(B U - B ld )(*j - *<) + ($,• - $0(^2, - C 2j )Y 

and considering the notations in (8) with l^andCi 
for 1 < i < r in (37)- (39), we can verify that the 
matrix inequality in (44) can be written as 



(45) 



In the same way, it is easy to check as the 
inequality(33) can be written as 



diag(T ^ , /, if, /, /, I) x (12) x diag(T \ , /, X2, /,/,/, ) <0 2. Note that a pole placement for each linear mod- 



(46) 



els of (TS) model in LMI region does not allow 
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specifying exactly the dynamics of the nonlinear 
model in closed loop. However, it is still a very 
interesting way of tuning the close loop dynam- 
ics. 



4 Numerical example 



values. The uncertain matrices are the following 



/ 0.06 0 0 \ 
H = 0 0.06 0 

\ 0 0 0.06/ 



(-G aT s 0 \ 

E 1 =E 2 =\ r 0 -M 2 T s x 10 -3 , 

V 0 0 -b J 

Ei i = E u = (0.1 0.1 0.1) , E 2l =E 22 = 0 



In this section, a numerical example is presented to 
illustrate the effectiveness of the proposed controller 
design technique. Consider the following discrete-time 
Lorenz system given in [7] with a sampling time T s = 
0.002s 



x(k + 1) = y>(xi(fc)) (Ai + AAi)x(k) 

i= 1 ^ 

+ (B, + A B u )u(k) + (Bj, + AB 2 , )w(k) J 

*(*) = E hi(xi(k))^C\ i x(k) + Di i u(fc)| 
y(k) = C 2 x(k ) + D 2 w(k) 

(51) 



where x( k) = [x±(k) xj ( k ) xJ(k)] T 
Similar to [7], we assume that xi (k) £ [Mi,M 2 ], 
and select the membership functions as hi = 
; h 2 = 1 — hi with Mi = -20 and M 2 = 30 
The disturbance input signal w(k) is given by w(k) = 
0.5sin(27rkTs). 

The system matrices are 



/l - aT s (jT s 0 \ 

Ai=\ rT s 1 - T s -MiT a , 

\ 0 MiT s 1 -bT s ) 

fl- vT s aT s 0 \ 

A 2 =\ rT s 1 - T s -M 2 T s , 

\ 0 M 2 T s 1-bTj 



B i 






C 2 




( o*i 

B 2l = -0.1 

V 0.1 
0 0 \ 
0.005 0 

0 0.005/ 





The nominal value of (cr, r, b) = (10,28, |) for chaos 
to emerge. We assume that all system parameters are 
uncertain but bounded within 30% of their nominal 



In this example, we choose the level 7 = 1 and 
we would place all its closed-loop poles within LMI 
disk region of radius ro =0.6 and center Co(— go, 0) = 
(-0.4,0). 

In order to design a fuzzy Hoo output feedback con- 
troller for the (TS) model, we first choose en = 
1.5, e 12 = 1.2, £22 = 1.2 and s 0ll = 1.3, e 0l2 = 

1.1, £022 1*1 

Then, using the Matlab LMI Control Toolbox with 
starting point xo = [10— 10— 10 ] and the theorem 
(3.2), we obtain 

f 55.396 21.341 -26.53l\ 

X = 21.341 25.134 -3.052 , 

\-26.531 -3.052 28.171 / 

/ 126.34 -43.164 4.6949 \ 

Y = -43.164 123.93 -0.35455 , 

\ 4.6949 -0.35455 120.36 / 

f 0.065816 -0.070841 0.31372 \ 

M = -0.091508 0.35792 0.072845 , 

y— 0.085573 0.15074 0.92246 / 

/ 0.097458 -0.16722 0.28155\ 

A 2 = -0.06961 0.34856 0.24871 , 

\— 0.073245 0.10855 0.90052/ 

/ 0.50427 -0.29683 \ 

B x = 0.17445 0.10098 

\0.048336 -0.071547/ 

/ 0.50256 -0.29948 \ 

B 2 = 0.17522 0.10008 

\0.04769 -0.072213/ 

Ci = (-0.29909 -0.15427 0.45631) , 

C 2 = (-0.29872 -0.15202 0.44942) , 

Using SVD algorithm, we have 

/— 73.713 -6.794 16.899 \ 

S= -17.48 -43.027-15.033 , 

\ 48.112 -26.042 20.428 / 

/ 76.39 11.157 -14.525\ 

W= 2.6467 44.994 14.07 
\— 47.027 20.656 -22.802 / 

Figure 1 shows the state responses and the ratio of 
the regulated output energy to the disturbance energy 

iYlk=o zT (^) z (^)/^k=o wT (^) w (^)) f° r chaotic 
system. We can note that the ratio tends to a con- 
stant value which is about 0.06488 which is less the 
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tion. IEEE transaction on automatic control, vol. 
42, no 7, July 1997. 

[3] C.S. Tseng, B.S.Chen and H.J.Uang. Fuzzy track- 
ing control design for nonlinear dynamic systems 
via T-S fuzzy model. IEEE transaction on auto- 
matic control, vol. 9, no 3, June 2001. 

[4] F.Liu, W.Xu and J.Chu. Robust fuzzy LQ-Hqo 
control for uncertain nonlinear Systems with de- 
lay 15th Triennial world congress, Barcelona- 
Spain, IFAC 2002 . 



Figure 1 : Control results of the chaotic 
Lorenz system. 







[5] H.D.Tuan, P.Apkarian and Y. Yamamoto. Pa- 
rameterized linear matrix inequality techniques 
in fuzzy contol system design. IEEE transaction 
on fuzzy systems, vol. 9, no 2 , April 2001. 

[ 6 ] H. Duong, P.Apkarian, T.Narikiyo and 
M.Kanota. New fuzzy control model and 
dynamic output feedback parallel distributed 
compensation. IEEE transaction on fuzzy 
systems, vol. 12, no 1, February 2004. 

[7] H.j.Lee, J.B.Park and G.Chen. Robust fuzzy con- 
trol of nonlinear systems with parametric uncer- 
tainties. IEEE transaction on fuzzy systems, vol. 
9, no 2 , April 2001. 



Figure 2 : Closed-loop poles of each sub- 
system . 



prescribed value 7 . 

Figures 2 shows that the fuzzy controller place all its 
closed-loop poles within prespecified LMI disk region. 



5 Conclusion 

In this paper, we have presented a stabilizing dy- 
namic output feedback controller which is devoted to 
uncertain nonlinear systems. The (TS) fuzzy model 
is employed to describe a nonlinear system with un- 
certainties. The Hqo performance and quadratic D- 
stability have been considered and a sufficient con- 
ditions, which are formulated in terms of LMI, is 
derived. Simulation results have clearly shown that 
the desired specifications for nonlinear systems can 
be achieved using the proposed method. 
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Abstract 

Under a deregulated environment, electricity consumers 
and suppliers generally establish various bilateral power 
transactions/contracts. The transmission company 
normally honors and executes these bilateral contracts as 
far as the system design and operating conditions permit. 
This paper describes studies of optimal the bilateral 
contracts by using Line Flow Factors (LFFs). A 
generalized linear programming formulation is then 
proposed to solve system operation under a deregulated 
environment subjected to the steady- state security 
constraints (e.g. generation and line flow limits). 
Examples are presented to illustrate how to use this 
formulation to minimize the cost of any bilateral contract 
to comply with the security requirements. It was 
concluded that the proposed methodology will be an 
effective tool to study the intricate relationships between 
the bilateral contracts and system security. 

Keywords: generation, load distribution factors and real 
power transaction. 

1. Introduction 

The open access transmission regime is spearheading the 
rapid disintegration of the well-entrenched vertically 
integrated structure of the electric power industry. The 
entry of a large number of new players and the 
unbundling of electricity services has pushed the industry 
toward the widespread use of transactions to meet 
customer demands. The driving forces of deregulation are 
aiming to establish a more competitive market in order to 
achieve lower rates for the consumers and higher efficiency 
for the suppliers. Traditional power companies are therefore 
gradually divested into independent business entities by 
unbundling and privatization of their generation, 

transmission and distribution functions [1,2]. 

Power suppliers, both the conventional utilities as well as 
the Independent Power Producers (IPP), are actively 
competing with one another for customers. The 
consumers can therefore establish various service 
contracts with any supplier in order to obtain the lowest 



rate and most desirable service. Bilateral contracts 
specifying the amount of power, the time and duration of 
the service and the associated rate and possible 
compensation are negotiated and agreed upon between 
the suppliers and consumers. 

The next step of course is to deliver the power from the 
suppliers to their respective consumers. Power 
transmission under a deregulated environment is 
generally handled by an independent entity whose 
network is open to all users. Although there are different 
rules governing the role and responsibility of the 
transmission company [2], a basic requirement is to serve 
the needs of all users in the network as much as the 
system design and operating conditions permit. In the 
competitive electrical power market, large-scale 
customers can select Power Producers and Suppliers 
(PPSs), and decide the amount of electrical power to 
purchase. This paper develops optimal bilateral 
transactions/contracts for customers for minimum cost 
taking into line limits. 

The line flow constraints are accounted for as a maiden 
attempt in determination of optimal bilateral transactions 
by expressing line flows in terms of active power 
contracts using line flow factors and are elegantly 
evaluated from existing load flow information. The 
proposed approach for optimal bilateral transactions with 
constraints on active power generations and line flows is 
tested on 4-bus sample system and 24-bus real life Indian 
system. The results obtained show great promise for 
practical application of the proposed algorithm on a real- 
time basis. 

The paper is organized as follows: Section (2) focuses on 
bilateral transaction model; section (3) presents the 
development of line flow factors, section (4) presents the 
evaluation of line flow factors and section (5) emphasizes 
the application of linear programming to obtain the 
optimal bilateral transactions. 
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2. Bilateral Contract/Transaction Model 

The electric power transactions are classified into pool 
transactions and bilateral contracts as well as security 
transactions. A certain transaction level of, for example, 
bilateral contracts may affect pool transactions by 
causing transmission congestion etc. Under the 
assumption that both the bilateral and market power 
transacted through pool (or power) market are coexisted 
and the generated power is supplied to customers through 
the transmission network management of power 
transmission company and the system operation of 
System Operator (SO). This paper proposes a method for 
determination of optimal bilateral contracts using linear 
programming with or without coexistence of pool market 
in a system for a given load condition. 

Bilateral transactions are contracts between power 
sellers (say GENCOS) and buyers (say DISCOS or large 
customers). This type of power trading would therefore 
entail injection of bulk power into the transmission 
system by the power producer, and withdrawal of an 
equal amount of power, from the network by the 
customer. This scenario is simulated as the output of a 
generator, at a generator bus, and a corresponding load, at 
a load bus. Following assumptions are resorted to: 

■ there are multiple candidate power suppliers, 
being included in the generation set; 

■ there are multiple candidate power customers, 
being included in the load set; 

■ the reactive power of load is compensated 
locally, only real power transaction is supplied 
by the transaction; 

The system model can be shown as in the following 
graph: 




< 






: Bilateral 
Contracts 



Figure 1 : The model for the bilateral transaction 

The goal of this paper is to find the optimal supplier 
customer combination for a particular transaction, such 
that the cost of transaction will have minimum value. The 
algorithm achieves these optimal bilateral transactions by 
considering all the transactions at a time. As the focus of 
the algorithm in this paper is on the determination of 
optimal bilateral power contracts, the actual flow path is 
of less importance. Therefore, a loss formula is preferred 
to the accumulation of all line losses. 

This paper focuses on the power dispatch problem under 
open access and is confined to the dispatch of active 
power, keeping bus voltages and reactive power constant 



when appropriate, The distributed slack bus concept 
where the amount of real power imbalance in the system 
can be distributed among transactions based on 
participation factors [3] resembles actual operation of 
power system and is easy to include but this is not done 
here to keep the presentation simple. 

The market design in the paper envisages participants 
who specify hourly quantities of demand, locations of 
generators supplying individual and group contracts, and 
willingness to pay for transmission services. At this stage 
the model does not take into account the co-existence of 
these bilateral and multilateral contracts with pool type 
dynamic supplies and demands based on bids and market 
clearing prices. This more complex issue of reconciling 
pool and bilateral contracts in an open access system is 
discussed in outline in [4-6]. Charging for Vars and other 
ancillary services is also an extension that stands outside 
the scope of this paper. 

3. Development of line flow factors 

To illustrate the proposed approach, a sample four-bus 
network shown in Fig.2 is considered. The system has 

two generators/suppliers, U and P '2 at buses 1, 2 and 

two customers/loads, Pd 3 and Pd 4 at buses 3, 4 
respectively. The system data is given in Appendix. 

Fet 

P gi : Supplier at bus i, 

P dj : Customer at bus j 

P d j : Power contract of customer at bus j from 
supplier at bus i 

P l _ k : Actual power flow in line 1-k due to all 
customers and suppliers 

P® k : Power flow in line 1-k due to supplier at bus i 

Pji k : Power flow in line 1-k with supplier at bus i and 
customer at bus j 

Pj_ k : Fine Flow Factor of a line 1-k when the 
supplier at bus i and customer at bus j 







Figurez. sample <4- dus systeinror the sample neiworK, h is assumed 
that the loads at buses 4 and 5 are supplied by both the suppliers. 
Therefore, we can write 
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p — p&t _i_ pg 2 
1 d 3 1 d3 1 d3 

p <,- p S+p£ (i) 

Where 

PH , P d \ = Power contracts of customers at buses 3 
and 4 from supplier at bust respectively. 

P d3 , P/ 4 2 = Power contracts of customers at buses 3 

and 4 from supplier at bus 2 respectively. 

Now the actual power flow in the line 1-2 can be 
expressed as 

P t -,= P&+P& ® 

Where 

P g _ 2 5 P g _ 2 = Power flow in the line 1-2 due to power 
contracts from supplier 1 and supplier2 respectively. 

The power flow P g \ can be expressed again as 



pg 1 _ p 13 I p 14 
1 1-2 1 1-2 ~ r 1 1-2 



(3) 



where 

P\\ > P\-i = P° wer flow in the line 1 -2 due to power 
contracts of customers at bus 3 and bus 4 from supplier at 
bus 1 respectively. 

Re-writing the above equation as 

P& = ( P t-2 / PS ) PS +( p t-2 1 p S ) PS 

p & HLF") Pg+( LF“) p S (4a) 

Where 

LF™ 2 = Line Flow Factor of line 1-2 or fraction of power 

flow in line 1-2 when supplier 1 contributes power P d3 
to customer at bus3. 

LF { [ 4 2 = Line Flow Factor of line 1-2 or fraction of power 
flow in line 1-2 when supplier 1 contributes power P d ^ 
to customer4 

Similarly the power flow P g _ 2 can be expressed as 

p il -( p “i p £ > p £ +( p,” / pS ) p £ 
p il -< £- 2 ) p £ +( £- 2 ) p £ + ( Si * ) pS m 

where 

L^l 2 = Line Flow Factor of line 1-2 or fraction of power 
flow in line 1-2 when the supplier at bus 2 contributes 
power P/ 4 2 to customer at bus 4. 

L\ 5 _ 2 = Line Flow Factor of line 1-2 or fraction of power 
flow in line 1-2 when the supplier at bus 2 contributes 
power P/ 5 2 to customer at bus 5. 

Similarly the power flow P g _ 2 can be expressed as 

p ll HP?,! p £) p £+( p£j p£ ) p £ 



p ll=( p F,%) p SHLF*) p S 

where 



(4b) 



LF^ 3 2 = Line Flow Factor of line 1-2 or fraction of power 
flow in line 1-2 when supplier2 contributes power Pf 3 2 to 
customer3. 

LF 2 2 = Line Flow Factor of line 1-2 or fraction of power 
flow in line 1-2 when the supplier2 contributes 
power P/ 4 2 to customer 4. 

From the equations 4(a) and 4(b) 

P — pg 1 _L pg 2 
1 1-2 1 1-2 1 1-2 

= (LF^PSMLF^pj; 

HLFS)Pg + aF? 2 )P£ (5) 

Similarly the actual power flow in the remaining lines 
have been expressed and represented in matrix form as 



> 1-2 " 




L u 

F-2 


L 14 

-^ 1-2 


Z 23 

^ 1-2 


T} 4 

M -2 




' P?' ” 


P\-3 




t\3 

M-3 


r } 4 

M-3 


r23 

M-3 


[} 4 

M-3 




1 d 3 
pg 1 


P\-A 


= 


z 13 

P-4 


t} 4 

-M-4 


L 23 

M-4 


I 24 

M-4 




r d4 

p £ 

pg 2 

_ r d4 _ 


p u 




t } 3 

^2-4 


L u 

^ 2-4 


L 23 

^ 2-4 


t} 4 

•^2-4 




>3-4 _ 




L u 

P-4 


t } 4 

^3-4 


L 23 

^3-4 


t} 4 

-^ 3-4 





In general, the relationship between line flow and the 
transaction power at a bus via LFFs can be expressed as 
follows: 

P,- k =^i-MU (7) 



4. Evaluation of line flow factors 

The proposed approach for calculating Line Flow Factors 
starts from the power flows, as given by a load flow 
solution, with the assumption that the network consists of 
one supplier and one customer at a time, as illustrated in 
Figure .2. A dominion is obtained for each supplier and 
each customer and is a directed graph consisting of one 
supplier and one customer. The Line Flow Factors for the 
customer at a bus j and supplier at bus i are calculated by 

dividing the real power flow in a line P } _ k by the load 
P dj at the bus j i.e., Fj_ k = P \_ k / P dj .The value of P dj is 

any arbitrarily selected value and need not be the actual 
value. If the given system is a lossless system then 

P dj = P gi otherwise P gi > P dj ■ 

The coefficient F\_ k , which is the sensitivity of the 

flow on the line 1-k with respect to the supplier at bus -i 
and customer at bus -k, implies the fraction of power 
flow in the transmission line 1-k when a unit power 
transaction is at bus-j . 

For example, if we assume that the given system is 
lossless, and that 100MW is the only customer at bus3 
(P d3 ) and the supplier at busl, then the line flows are 

obtained as given by a load flow solution and the line 
flow factor of line 1-2 for this dominion is 
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L\\ =16.35/100 =0.1 635. Here the line flow 

P x _ 2 ( = Pi-k ) i n the line 1-2 is 16.35MW and P d3 ( : =P dj ) 

is 100MW. The subscript 1-2 indicates the line connected 
between buses 1 and 2 and the superscript 1-3 indicates 
that supplier is at bus 1 and customer is at bus 3. 
Similarly the line flow factors for other dominions are 
calculated and are given in Table I. 

Removal of 




A, Ki 



Figure 3: A domain 



Table 1: Line Flow Factors of Each Dominion 







Supplier at busl 


Supplier at bus2 


Lines 


LUFs 


(i=l) 


(1=2) 


1-k 




Customer at bus 


Customer at bus 


n 1 


3 


4 


3 


4 




^ i~k 


0 = 3 ) 


0 = 4 ) 


0 = 3 ) 


0 = 4 ) 


1-2 


A-2 


0.1635 


0.3000 


-0.3052 


-0.1687 


1-3 


A-3 


0.6180 


0.3007 


0.3906 


0.0725 


1-4 


A-4 


0.2177 


0.3993 


-0.0854 


0.0962 


2-4 


^2-4 


0.1635 


0.3000 


0.6949 


0.8313 


3-4 


^3-4 


-0.3811 


0.3007 


-0.6093 


0.0725 



The negative sign indicates that, the real power flow is 
opposite to the convention used. Once the line flow 
factors are obtained, these factors remain constant, and 
can be used for calculating line flows without running 
load flows whenever load/generation at a bus is increased 
or decreased. 

5. Optimal Bilateral Transactions/Contracts 
(OBT) Problem Formulation 

The optimal bilateral transactions problem is considered 
as a general minimization problem with constraints, and 
can be written in the following way: 



Min Z = f T x (8) 

Such that A x<b 

A eq X ^ b e q 

lb < x <ub (9) 

where f , x , b , b eq , lb , and ub are vectors and 

A and A pn are matrices. 

eq 

For the four- bus test system the cost coefficients 

/ x 

vector — and — is given by: 



"15" 




Ug l 
r d 3 


15 




pfi 




, x = 


17 




p£ 


_ 17 _ 




1 

K> 

1 



The elements of the vector f represents marginal cost 

coefficients, and the elements of the vector x represents 

the bilateral contracts between each supplier and 
customer. 

The objective function of our linear program is the sum 
of marginal cost of real power suppliers at buses 1, and 2. 
We denote the objective function of the linear program as 
Z: 

Z=15P/ 3 1 + 15* + 11 pf, + 17P/ 4 2 (10) 

Equality Constraints: 

Denoting the total load to be supplied as LOAD, the 
original equality constraint is that the sum of the bilateral 
power contracts must equal the total load to be supplied, 
i.e., 

P// + P/ 4 ‘ + P/ 3 2 + P/ 4 2 = LOAD =500MW (11) 

Also the power supplied to the customer from a supplier 
should be equal to the load power i.e. 

P£+P£=P*,= 300MW, (12) 

Pja + Pf a = P</4 = 200MW (13) 

Inequality Constraints: 

The original inequality constraints are 

100<P gl <500, 50 < P g2 <150, 

However, in order to utilize 4 variables, the above 
constraints on the two-generation levels are converted to 
constraints on the 4 new variables. 
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6. Real-Life Indian System studies 

The proposed approach is applied for a practical Indian 
system of 24-bus equivalent EHV power system and is 
shown in Fig.4. The system has 4 generator buses and 20 
other buses. The load is represented at 220 kV sides of 
400 kV / 220 kV nodes at 8 numbers of buses. The 
system peak load is 2620MW; 980MVAR.lt is assumed 
that all the generators are having thermal units. . 



The equations (11)-(13) form the additional constraints 



represented in the matrix A eq 

Lower Bounds : It should be recognized that all of the 
new variables must be non-negative. This comes from the 
fact that we cannot have a negative amount of generation 
increment. But we can have any amount of generation 
increment between zero and the upper bound. So the 
lower bound on all 4 variables is zero. 

Upper bounds : The upper bounds on the variables can be 
taken as the maximum load at that particular load bus. 
For the four-bus system it is taken as 200 MW. 

The results of the linear programming are given by 



1 

a ? 

1 




' 200 . 0 ' 


ps l 
r d 4 




200.0 


ps 2 
r d 3 




100.0 


1 

K> 

1 




°. 0 _ 



The summary of the results obtained using the LFF 
method are given in the Tables 6. 7. and 6.8. 



Table 2: Real Power Transactions & Cost 



Real Power 
Contracts 


Transaction 

(MW) 


Cost of Transaction 
($/MWhr) 


Pfi (MW) 


200 


3000 


Pf 2 (MW) 


100 


1700 


>4 (MW) 


200 


3000 


Pd 4 (MW) 


0 


0 


Total Cost ($/MWhr) 


7700 



o« 




Table 4: Generating unit’s data 



Generator 
Bus No. 


No. of 
Units 


p 

gmax 


p 

gmrn 


1 


10 


200 


40 


2 


2 


200 


40 


3 


4 


200 


40 


4 


4 


200 


40 



Table 4: Cost characteristics of generating units 



Generator 
Bus No. 


a -MW 2 


b -MW 


c 


1 


0.8 


400 


8000 


2 


0.8 


400 


8000 


3 


0.8 


400 


8000 


4 


0.8 


400 


8000 



Table 3: Real Power Generation Levels & Cost 





Generation 

levels 

(MW) 


Cost of 
Generation 
($/MWhr) 


G1 (MW) 


400 


6000 


G2 (MW) 


100 


1700 


Total Cost ($MWhr) 


7700 



7. Conclusions 

In this paper, a new set of line flow factors are 
developed. In the developed framework, one can directly 
evaluate the flow on a line and optimal bilateral real 
power contracts. The method is conceptually simple. The 
linear programming concept is used to determine the 
optimal bilateral contracts, which minimizes the marginal 
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cost. For the sample 4-bus network and real -life 24-bus 
Indian systems, the method has been demonstrated with 
the encouraging results. 



Table 6: Real Power contracts 



Load 

Bus 

No. 


Real power contracts from each load each bus 
(MW) 


G1 


G2 


G3 


G4 


5 


187.2724 


46.8410 


163.0704 


36.5883 


6 


29.1313 


31.7595 


55.2295 


184.8875 


7 


42.3607 


34.1981 


80.3218 


186.2895 


8 


28.6762 


52.0084 


54.3807 


66.9826 


9 


24.9035 


18.2312 


47.2200 


39.6016 


10 


13.2460 


10.4365 


25.1163 


16.5207 


13 


73.0557 


106.3393 


138.5519 


180.5350 


15 


531.3374 


31.6870 


110.3275 


24.7549 



Table 7: Generation levels (MW) and cost (Rs./hr) 



Generator 


Method 


Cost of generation 
Rs./hr 


G1 


929.95 


521160 


G2 


331.53 


044280 


G3 


674.23 


115820 


G4 


736.16 


167130 


P. Loss (MW) 


51.87 


Total Cost: 848390 
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Appendix 



Table A: The four- bus system data 



Transmission line data 


Generation and Load data 


From 


To 


R 


X 


Bus 


pg 


Pd 


Bus 


Bus 


(CS) 


(Q) 


No. 


(MW) 


(MW) 


1 


2 


12.75 


97 








1 


3 


6 


69.5 


1 


400 


0 


1 


4 


11.7 


96 


2 


100 


0 


2 


4 


3.5 


30.8 


3 


0 


300 


3 


4 


5.75 


58 


4 


0 


200 
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Abstract 

The proportional-integral-derivative (PID) controllers 
were the most popular controllers of this century because 
of their remarkable effectiveness, simplicity of 
implementation and broad applicability. However, PID 
controllers are poorly tuned in practice with most of the 
tuning done manually which is difficult and time 
consuming. The computational intelligence has purposed 
genetic algorithms (GA) and particle swarm optimization 
(PSO) as opened paths to a new generation of advanced 
process control. These advanced techniques to design 
industrial control systems are, in general, dependent on 
achieving optimum performance with the controller 
when facing with various types of disturbance that are 
unknown in most practical applications. This paper 
presents a comparison study of using two algorithms for 
the tuning of PID-controllers for processes which 
represents a subsystem of complex industrial processes, 
known to be non-linear and time variant. Simulation 
results showed that the PID control tuned by PSO 
provides an adequate closed loop dynamic for the Ball 
and Hoop system experiment in wide range operations. 

Keywords: Particle swarm optimisation; genetic 

algorithms, intelligent control; PID control 

1. Introduction 

Currently, the most of industrial processes are controlled 
by the well-established PID control. The standard PID 
control configuration is as shown in Figure 1. These 
parameters K p , K { , K d are chosen to meet prescribed 

performance criteria, classically specified in terms of rise 
and settling times, overshoot, and steady state error, 
following a step change in the demand signal. The 
popularity of PID control can be attributed to its 
simplicity (in terms of design and from the point of view 
of parameter tuning) and to its good performance in a 
wide range of operating conditions. Moreover, the 
controllers designed with the aid of modern control 
techniques are usually of high order, difficult to 
implement, and virtually impossible to re-tune on line. 
Also, the design procedures of PID controllers are 



simple, although, the tuning parameters are the most 
challenging ones in a PID controller, as optimal tuning 
parameters are difficult to find. Therefore, researchers 
introduced advanced controllers such as PID auto-tuning 
or self-tuning, model-based adaptive control, model 
predictive control, fuzzy control, artificial neural 
networks, optimal control, expert control, robust control, 
etc. Although fuzzy control has great potential for 
solving complex control problems, it has, at least, these 
limitations 

1. Its design procedure is complicated and requires a 
great deal of specialty. 

2. As a rule-based system, criteria for control system 
stability are complicated and usually conservative. 
A bad-rule or an out-dated rule may cause the fuzzy 
controlled system to become unstable. 

3. If the process changes, it is necessary to modify the 
rules, which is time-consuming. 

The other methods including optimal control, expert 
control, neural networks, etc., have been developed to 
address the issues of complexity. However, most require 
some kind of process model and high-level expertise. It 
is difficult to develop a general-propose, user-friendly, 
smart control system based on these control methods. 




Figure 1: PID Controlled System. 



There are a variety of PID controllers tuning methods 
have been developed such as Ziegler-Nichols rules, 
Cohen-Coon rule and so on [1, 2, 11]. These methods are 
applied directly since they provide simple tuning rules to 
determine the PID parameters. However, since they rely 
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on a minimum amount of dynamic information by 
making a certain assumption about nature of the 
controlled process, such as linearity, weak interactions 
within the process, absence of noise, etc. The realized 
closed loop response is less than optimum since the real 
world processes are non-linear and very complex. 
Recently, the computational intelligence has purposed 
genetic algorithms (GA) [9, 18] and particle swarm 
optimization (PSO) [15, 16] as opened paths to a new 
generation of advanced process control. These advanced 
techniques to design industrial control systems are, in 
general, dependent on achieving optimum performance 
with the controller when facing with various types of 
disturbance that are unknown in most practical 
applications. Usually such an optimum system 
performance is defined by a continuous function which 
is obtained by minimizing the sum of the square of the 
error between the actual and desired output. The GA and 
PSO are non-gradient based optimization and search 
algorithm and belongs to probabilistic search algorithms, 
which have attracted much attention from the research 
community. These algorithms generally mimic some 
natural phenomena, for example GA model the evolution 
of a species, based on Darwin’s principle of survival of 
the fittest and competition to produce better adapted 
generations in their problem solution space while the 
PSO in contrast relies on cooperation rather than 
competition and good solutions in the problem set are 
shared with their less-fit brethren so that the entire 
population improves and poorly performing members are 
not killed off as in GA. Although these probabilistic 
search algorithms generally require many more function 
evaluations to find an optimal solution, as compared to 
gradient-based algorithms (i.e., simple gradient descent, 
conjugate gradient descent, etc.), they are promising 
approaches due to their effectiveness in searching very 
large spaces and the ability to perform global search for 
best forecasting model. Moreover, they are generally 
easy to program, can efficiently make use of large 
numbers of processors, do not require continuity in the 
problem definition, and generally are better suited for 
finding a global, or near global, solution. Genetic 
algorithms have been extensively applied to the off- 
line/on-line design of controllers. GA was first 
investigated as an alternative means of tuning PID 
controllers. Oliveira et al. [19] used a standard GA to 
determine initial estimates for the values of PID 
parameters. Genetic algorithms are tested for PID 
controller for nonlinear process and showed robustness 
and efficiency [22, 12, 10]. Lennon and Passino [17] 
developed a more complex optimization known as a 
general genetic adaptive controller (GGAC). Other 
several applications of GA appeared in [4, 7, 21]. The 
PSO have also been used for tuning PID controller [8, 
20; 6, 23] and have proved successful results. The 
primary contribution of this paper is to evaluate the two 
algorithms in the tuning of PID-controllers in the ball 
and hoop system which represents a system of complex 
industrial processes, known to be non-linear and time 
variant. Such comparative analysis is very important for 
identifying both the advantages and their possible 
disadvantages. The rest of the paper is organized in the 
following manner. In Section 2, a problem formulation is 



presented. In sections 3 and 4, the particle swarm 
optimization algorithm and genetic algorithms which 
will be applied in this paper are described briefly. In 
section 5, we provide an overview of control problem 
which we will be used in this paper. Detailed 
experiments results are provided in section 6. Finally, 
Section 7 concludes the paper. 

2. Problem Formulation 

By formulating the PID controller design problem as an 
optimization and search problem, we frequently have a 
complicated function to be optimized. This function can 
be represented as follows: 

f-S^R ( 1 ) 

Where S is the set of solutions and best choices are that 
for the function / is optimal. In the context of controller, 
a candidate system can be represented by a uniform 
parametric vector given by 

Si={g\,g2,---,g n }^ n ( 2 ) 

Where i stands for the i th possible candidate solution, n 
the number of parameters required by the solution, 
gj e 91 the 7 th parameter of the i th candidate solution 

withy e {l,...,«}, and 91 n the ^-dimensional real 
Euclidean space. We have a complicated objective 
function g(5) needs to be optimized. Where Q(S) 
represents the quality measurement for a solution Sf 
given VS/ g(s z - ) > 0 . The problem is to find the best 
solution (i.e., controller) S such that: 

Q{s)=MaxQ(s) (3) 

s 

3. Particle Swarm Optimization 

The particle swarm optimization algorithms are based on 
two socio-metric principles. Particles fly through the 
solution space and are influenced by both the best 
particle in the particle population and the best solution 
that a current particle has discovered so far. The best 
particle in the population is typically denoted by (gobal 
best), while the best position that has been visited by the 
current particle is donated by (local best). The (global 
best) individual conceptually connects all members of 
the population to one another. That is, each particle is 
influenced by the very best performance of any member 
in the entire population. The (local best) individual is 
conceptually seen as the ability for particles to remember 
past personal success. The particle swarm optimisation 
makes use of a velocity vector to update the current 
position of each particle in the swarm. The position of 
each particle is updated based on the social behaviour 
that a population of individuals adapts to its environment 
by returning to promising regions that were previously 
discovered [14]. Let the z'-th particle of the swarm is 
represented by the D-dimensional 

vector X/ = (*/i,*/2 >—>*/£)) and the best particle in the 
swarm, i.e. the particle with the smallest function value, 
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is denoted by the index g. The best previous position (the 
position giving the best function value) of the z'-th 
particle is recorded and represented as 
Pf = {pibPi2>-,PiD)> an d the position change (velocity) 
of the 7-th particle is Vf = (v/i,v/2 v-jV/d) • The particles 
are manipulated according to the equations 



v id + C 1 -n ip id - x id ) 


( 4 ) 


+ °2 -/'2 \Pgd ~ x id ) 


x id ~ x id + v id 


( 5 ) 



where d = 1, 2, . . . , D\ i = 1, 2, . . . , N and N is the size 
of population; w is the inertia weight; c\ and c 2 are two 
positive constants; r x and r 2 are two random values in the 
range [0, 1]. The first equation is used to calculate z'-th 
particle’s new velocity by taking into consideration three 
terms: the particle’s previous velocity, the distance 
between the particle’s best previous and current position, 
and, finally, the distance between swarm’s best 
experience (the position of the best particle in the 
swarm) and z-the particle’s current position. Then, 
following the second equation, the z'-the particle flies 
toward a new position. In general, the performance of 
each particle is measured according to a predefined 
fitness function, which is problem-dependent. The role 
of the inertia weight w is considered very important in 
PSO convergence behaviour. The inertia weight is 
employed to control the impact of the previous history of 
velocities on the current velocity. In this way, the 
parameter w regulates the trade-off between the global 
(wide-ranging) and local (nearby) exploration abilities 
of the swarm. A large inertia weight facilitates global 
exploration (searching new areas); while a small one 
tends to facilitate local exploration, i.e. fine-tuning the 
current search area. A suitable value for the inertia 
weight w usually provides balance between global and 
local exploration abilities and consequently a reduction 
on the number of iterations required to locate the 
optimum solution. A general rule of thumb suggests that 
it is better to initially set the inertia to a large value, in 
order to make better global exploration of the search 
space, and gradually decrease it to get more refined 
solutions, thus a time decreasing inertia weight value is 
used. The main steps of the PSO algorithm are shown in 
Table 1, where there are three main steps 1) initialize a 
population of particles (position and velocities); 2) 
updating velocities; 3) updating positions. 

Table 4. The main steps for PSO algorithm. 

Function PSO(Particles, Fitness-FN) returns a controller 

inputs: Particles, a set of PID-controllers 

Fitness-Fn, a function that controller response according Eq. (7). 

repeat 

Loop for i from I to Size (Particles) do 

Find the personal best S global best position of each Particle. 
Update each Particle according to Eq. (4) 5 Eq. (5). 
until some controller is fit enough, or enough time has elapsed 
return the best controller in Particles, according to Fitness-Fn 



4. Genetic Algorithm 

The GA is an optimization routine based on the 
principles of Darwinian Theory and natural genetics. It 
has primarily been utilized as an off-line technique for 
performing a directed search for the optimal solution to a 
problem. The GA performs a parallel, directed, random 
search for the fittest element of a population within a 
search space. The population simply consists of strings 
of numbers, called chromosomes that hold possible 
solutions of a problem. The members of a population are 
manipulated cyclically through three primary genetic 
operators called selection , crossover , mutation, and 
replacement, to produce a new generation (a new 
population) that tends to have higher overall fitness 
evaluation. By creating successive generations which 
continue to evolve, the GA will tend to search for a 
global optimal solution. The key to the search is the 
fitness evaluation. The fitness of each of the members of 
the population is calculated using a fitness function that 
characterizes how well each particular member solves 
the given problem. Parents for the next generation are 
selected based on the fitness value of the strings. That is, 
strings that have a higher fitness value are more likely to 
be selected as parents, and, thus, are more likely to 
survive to the next generation. The main steps of the GA 
algorithm are shown in Table 2. 

Table 4. The main steps for GA algorithm. 

Function GA(Pnpulatinn, Fitness-FN) returns a controller 

inputs: Population, a set of PID-controllers 

Fitness-Fn, a function that controller response according Eq. (7). 

repeat 

New-Population <- empty set 
Loop for i from I to Size (Population) do 

Father <- Random-Selection (Population, Fitness-Fn) 

Mother <- Random-Selection (Population, Fitness-Fn) 

Child <- Crossover (Father, Mother) 
if (small random probability) then 
Child <- Mutate (child) 
add Child to New-Population 
Population <- New-Population 

until some controller is fit enough, or enough time has elapsed 
return the best controller in Population, according to Fitness-Fn 

5. Plant System 

The Ball and Hoop system [3] illustrates the dynamics of 
a steel ball that is free to roll on the inside of a rotating 
circular hoop. There is a groove on the inside edge of the 
hoop so that a steel ball can roll freely inside the hoop. 
This introduces the complexity of the rolling radius of 
the ball being different to the actual radius of the ball as 
illustrated in Figure 2, where the angle, 6 , is the hoop 
angular position. The position of the ball is given by: 

1. y, the position of the ball on the hoop periphery 
with respect to a datum point. 

2. y / , the slosh angle which measures the deviation of 
the ball from its rest position. 
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A fourth order system was selected as the Ball and Hoop 
system is of order four with the following transfer 
function [10]. 



4v)= 



i 

j 4 +6-s 3 +11-s 2 +6-j 



( 6 ) 



The ball and hoop apparatus is difficult to control 
optimally using a PID controller because the system 
parameters are constantly changing. The Ziegler-Nichols 
tuning method [24] using root-locus to evaluate the PID 
gains for the system, and the resulted gains are: 
K p =6,Ki =1.91, and Kj =4.74. 




6. PID Process Control Design 

The PID controller will be tuned offline using both of 
GA and PSO as shown in Figure 3. 



dimension. The velocity of each particle is randomized to 
a small value to provide initial random impetus to the 
swarm. The swarm size was limited to 20 particles. The 
must important factor is maximum velocity parameter 
which affect the convergence speed of the algorithm is 
set to 100.0. The c x and c 2 are 2.0 and 2.0 respectively. 

For GA, the GAOT is used [13] with population size of 
50, the search range [-100.0, 100.0] and with other 
default parameters. The two algorithms are runs of 100 
iterations. Figure 3 presents the step response of different 
PID controller using ISE performance index. As shown 
in Figure 4, the best results when used in conjunction 
with a particle swarm optimization. Table 3, describes 
the steady state characteristics of each of the controlled 
systems. 




Figure 4: The step response of different tuning algorithms. 




Figure 3: The structure of GA or PSO of PID tuning system. 



Where the performance index ISE is used to estimates 
the parameters of PID controller is given by: 

T 

I ISE = Je 2 (/)A (7) 

0 



Table 1 : Comparisons of Steady State Responses 





PSO 


GA 


ZN 


Rise Time 


1.2000 


1.2000 


2.1000 


% Overshoot 


24.345 


28.5908 


58.1561 


Settling Time 


10.0 


20.4000 


19.8000 



The P SO-tuned PID controllers outperform the GA- 
tuned and Ziegler-Nichols tuned controller in terms of 
settling time and overshoot. However the rise time of 
both PSO and GA are equal. Analyzing the transfer 
function of the ball and hoop system is stable as it has 
three poles located on the left hand side of the s-plane 
and one critically stable pole located at the origin. The 
presence of a critically stable pole will result in a more 
oscillatory open loop response. From the transfer 
function, it is clear that there will be an open-loop steady 
state gain of 1.84. This gain is shown in Figure 5 when 
the open loop response of the system is plotted against a 
PID controlled system using the PSO tuned controller 
which gave the best performance. 



The tuning of PSO-based PID controller' results will be 
compared to those obtained from the GA and traditional 
techniques. The PSO and GA tuning algorithms are 
conducted several pre-experiments to determine the 
parameters setting per algorithm that yields the best 
performance. For PSO, all swarm particles start at a 
random position in the range [-100.0, 100.0] for each 



Figure 6 compares the performance of the Ziegler- 
Nichols tuned designed controller with the PID tuned 
controller by the PSO and GA. The ISE performance 
criteria were improved by 35% for the GA and 52% for 



the PSO. 
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Figure 5: The open and PID closed loop response of the Ball and Hoop 
system tuned by PSO. 



Tunining PID Controller 




Time (sec) 

Figure 6: The Ball and Hoop System 



7. Control System Stability Criterion 

In control systems, which is a linear time-invariant 
system its stability is determined by its characteristic 
polynomial and its corresponding roots, which are 
eigenvalues or poles. A linear time-invariant system is 
asymptotically stable if and only if all the roots of its 
characteristic polynomial are in the open left half of the 
s-plane. The potential for closed-loop pole location is 
determined by the flexibility available to compute 
acceptable three-term controller coefficients to achieve 
arbitrarily chosen stable closed-loop pole positions. For a 
discrete system the poles of the overall closed loop 
equation must lie with the unit circle. Any set of PID 
parameters that gave poles that were outside the unit 
circle had there weighted error value set really high so 
that they would not be chosen for selection in the GA or 
PSO learning process. As the poles obtained for the 
overall closed loop transfer function for the different sets 
of PID parameters can be real and imaginary the absolute 
value of each pole is taken. If this absolute value is 
greater than one then the pole is outside the unit circle 
hence the overall system is unstable. When the system is 
unstable for a set of PID values the fitness evaluation 
gives this set of PID parameters a low fitness value. 



8. Discussion and Conclusions 

In this work, we have presented a comparison study of 
using PSO and GA algorithms for the tuning of PID- 
controllers for processes which represents a subsystem of 
complex industrial processes, known to be non-linear 
and time variant. Simulation results showed that the PID 
control tuned by PSO provides an adequate closed loop 
dynamic for the Ball and Hoop system experiment in 
wide range operations. This work can be extended by 
tuning the PID controller online so that the controller 
adaptively changing its setpoint based on the current 
situation. 
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ABSTRACT 

Synchronous static compensator (STATCOM) can be 
used to improve the dynamic performance of a power 
system. This article presents an online adaptive pole- 
shift method for stabilization of a single machine 
system which is then extended to a multi-machine 
system. An adaptive linear plant parameter model is 
used to derive the pole-shift control strategy. The 
controller performance has been tested for various 
disturbances. Responses for the single machine 
system were also compared with optimized PI control 
and robustness of the proposed adaptive algorithms 
was tested. From a number of simulation studies on 
the single machine as well as the multi-machine 
power system it was observed that the adaptive 
algorithm converges very quickly and also provides 
robust damping profiles to the under damped power 
system. 

Keywords 

Power system stabilizing control, STATCOM, on-line 
identification, adaptive control, pole - shifting control 



Nomenclature 



8 


Generator rotor angle 


CO 


Rotor speed 


Vo 


Base (synchronous) speed 


Pm 


Mechanical power input 


Pe 


Electrical power output 


M 


Inertia constant 


D 


Damping coefficient of generator 


e q , 


Quadrature (q) axis internal voltage 


Tdo 


Open circuit field time constant 


Efd 


Field voltage 


X d , X d 


Synchronous, transient direct (d) axis 




reactance 


Id 


d-component of armature current 


K a , T a 


Exciter gain, time constant 



V t Generator terminal voltage 

V L STATCOM bus voltage 

V b Network bus voltage 

Efdo, V to Nominal field, terminal voltage 
Vdoldc dc capacitor voltage, current of 

STATCOM 

C dc Capacitance of dc capacitor 

I Lo Steady ac STATCOM current 

m, \j/ Modulation index, phase of STATCOM 
voltage 

0 Parameter vector 

cp Measurement Vector 

X Forgetting Factor Bus admittance 

matrix of system excluding generator 
and STATCOM 

1. Introduction 

The static synchronous compensator (STATCOM) is 
a power electronic based synchronous voltage 
generator that generates a three-phase voltage from a 
DC capacitor. By controlling the magnitude of the 
STATCOM voltage, the reactive power exchange 
between the STATCOM and the transmission line and 
hence the amount of shunt compensation in the power 
system can be controlled [1]. In addition to reactive 
power exchange, properly controlled STATCOM can 
also provide damping to a power system [2, 3]. 

The modeling, operation and control fundamentals of 
the STATCOM have been extensively discussed in 
the literature [1, 4-7]. While most of the control 
designs are carried out with linearized models, 
nonlinear control strategies for STATCOM have also 
been reported recently [8]. STATCOM controls for 
stabilization have been attempted through complex 
Fyapunov procedures for simple power system 
models [8]. Applications of fuzzy logic and neural 
network based controls have also been reported [9, 
10]. Stabilizers based on conventional linear control 
theory with fixed parameters can be very well tuned 
to an operating condition and provide excellent 
damping under that condition. But they cannot 
provide effective control over a wide operating range 
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for systems that are nonlinear, time varying and 
subject to uncertainty. In order to yield satisfactory 
control performance, it is desirable to develop a 
controller which has the ability to adjust its 
parameters from on-line determination of system 
structure or model, according to the environment in 
which it works. Application of adaptive control 
methods to power system excitation control systems 
and static VAR systems has been reported [11, 12]; 
however, its application to FACTS devices has been 
very limited. 

This article investigates the stability enhancement 
problem of power systems installed with STATCOM 
through adaptive online control of converter voltage 
magnitude. The adaptive control strategy developed 
has been tested for its robustness over wide ranges of 
operation for a single machine as well as a multi- 
machine system. 

The organization of the article is as follows: Section 2 
gives the model of the power system including 
STATCOM; the theory of self-tuning adaptive 
regulator is included in section 3, and section 4 
presents the pole-shift control algorithm used in 
conjunction with adaptive control. Sections 5, 6, and 
7 devote on the simulation results for single machine 
case, while sections 8 and 9 give simulation results on 
the multi-machine system. At the end the conclusions 
of the study are presented in section 9. 

2. Single Machine System with 
STATCOM 

A single machine infinite bus system with a 
STATCOM installed at the mid-point of the 
transmission line is shown in figure 1 . The 
STATCOM consists of a step down transformer, a 
three phase GTO-based voltage source converter, and 
a DC capacitor. The STATCOM is modeled as a 
voltage sourced converter (VSC) behind a step down 
transformer. The dynamic equations of the generator- 
excitation system are [13], 




5 = co 0 Aco 

co — — — [P m -P e -DA co\ 

M 

e q =^-[E fd -e q -(x d -x d )I d ] (1) 

1 do' 

Efd =-T(E fd -E fd0 ) + ^(V t0 -V t ) 

A A a a 

The state variables [8 co e q E fd ] are the generator 
rotor angle, speed, internal voltage and field voltage 
respectively. The STATCOM capacitor voltage 
equation is, 

(2) 

Here, m and \| / are the magnitude and phase angle of 
the STATCOM ac side voltage. Combining (1) and 
(2), the dynamics of the single machine system is 
written as a 5 th order state model as, 
x = f(x,u) (3) 

Here, the control is considered to be the vector 
comprising of m and \p, respectively. 

3. Self-tuning Adaptive Regulator 

Self-tuning control is one form of adaptive control 
which has the ability of self-adjusting its control 
parameters according to system conditions. Figure 2 
shows the structure of an adaptive regulator. The self- 
tuning strategy is composed of two processes - the 
system identifier and the controller. The identifier 
determines the parameters of the mathematical model 
of the system from input-output measurement of the 
plant. The parameters of the identifier are updated at 
each sampling instant so that it can track the changes 
in the controlled plant. The control for the plant is 
then calculated based on this recursively updated 
system model. 




Figure 2 Block diagram of self-tuning controller 
The plant model is assumed to be of the form, 

A(z" 1 )y(t)=B(z" 1 )u(t)+e(t) (4) 




56 




Automatic Control and System Engineering Journal, ISSN: 1687-4811, Volume 7, Issue 2, ICGST, Delaware, USA, Nov. 2007 



where, y(t),u(t) and e(t) are system output, input and 
the white noise, respectively; z" 1 is the delay operator. 
The polynomial A and B are defined as, 

A(z l ) = 1+ aiz 1 + a 2 z 2 + a 3 z 3 + a 4 Z 4 + ... (5) 

B(z _1 ) = 1+ biz" 1 + b 2 z" 2 + b 3 z" 3 + b 4 z" 4 + (6) 

The vector of parameters 0(t)= [ai a 2 ,bi b 2 

] T are calculated recursively on-line through the 

recursive least square [11] technique using, 

6(t+l) = 0(t) +K(t) [y(t) - 0 T (t) ¥ (t)] (7) 

The measurement vector, modifying gain vector, 
and the covariance matrix, respectively are, 



v(t)= [-y(t-l) y(t-2) y(t-n a ) u(t-l) u(t-2). 

n b )] T 



m= 



h(t) + y/ T (t)P(t)y/(t) 



P(t + 1) 



( 8 ) 



1 

m 



[P(t)-K T (t)P(tV(t)] 



•u(t- 



T(z _1 )= A(z l )¥(z l ) + B(z l )G(z l ) (10) 

The pole-shifting algorithm makes T(z _1 ) take the 
form of A(z _1 ) but the pole locations are shifted by a 
factor a, i.e. 

ACz ^FCz 1 ) + B(z ')G(z ') = A(a z 1 ) (11) 

Expanding both sides of (11) and comparing the 
coefficients give, 



' 1 


0 


. 0 


bi 


0 


. 0 " 


" f, " 




a^a-l) 


a l 


1 


. 0 


b 2 


b l 


. 0 






a 2 (a 2 -1) 




a l 






b 2 


. 0 








a n 

n a 




. 1 


b n b 




• b b 


f„ f 






0 


a n 

n a 


. a 1 


0 


b n b 


. b 2 


go 




a na (« na -l) 




0 






0 








0 


0 


0 


. a n 

n a 


0 


0 


• V. 


>g. 




0 



The above is written in the form, 

M Z(a )=L(a) (12) 



X(t) is the forgetting factor; n a and n b denote the order 
of the polynomials A and B, respectively. The 
identified parameters in (7) can be considered as the 
weighted sum of the previously identified parameters 
and those derived from the present signals. 

4. The Control Strategy 

Using the parameters obtained from the real time 
parameter identification method, a self-tuning 
controller based on pole shift is computed on-line and 
fed to the plant. Under the pole shifting control 
strategy, the poles of the closed loop system are 
shifted radially towards the centre of the unit circle in 
the z-domain by a factor a, which is less than one. 
The procedure for deriving the pole-shifting 
algorithm [14] is given below. 

Assume that the feedback loop has the form, 

«(0 _ <j(z -1 ) (9) 

y(t) F{z~ x ) 

where, 

F(z _1 ) = 1+ f x z~ l + f 2 z" 2 + f 3 z" 3 + f 4 z" 4 + 

+fnfZ' nf 

G(z‘) = g 0 + giz' 1 + f 2 z' 2 + f 3 z' 3 + f 4 z' 4 + +f ng z' 

ng 

n f = n b -l, n g = n a -l 

From (4) and (9) the characteristic polynomial can be 
derived as, 



If parameters [{ai}, {b)}] are identified at every 
sampling period and pole-shift factor a is known, the 
control parameters Z=[{fi},{gi}}] solved from (12) 
when substituted in (7) will give, 

u(t, a) = X T (t)Z = X T (t)M‘ l L(a) ( 1 3) 

In the above, X(t)=[-u(t-l) -u(t-2) u(t-n f ) -y(t) -y(t- 

1) -y(t-2) y(t-n g )] 

The controller objective is to force the system output 
y(t) to follow the reference output y r (t). The objective 
function can then be expressed as, 

J = mm[y(t)-y r (t)] 2 (14) 

a 

In the above, y(t)=biu(t) + X T (3; (3=[-b 2 -b 3 

ai a 2 ]. 

If the variation of J with respect to a can be made 
smaller than a predetermined error bound 8i ? it can be 
shown that a minimum of J will occur when, 

Aa = JizM ( i5) 

£ 2 +~[fJ 3 +2b?f 2 2 ] 

where, s 2 is a small number chosen to avoid the 
singularity. In the above, 

. _ 8 J _ du d 2 u 

J l — " 7 T ? J 2 ~ T ? J 3 — ~Z T 
ou da da 

The partial derivatives are evaluated from (13) and 
(14), and updates of the control is obtained 
considering first few significant terms of the Taylor 
series expansion of u(t 9 a). The algorithm can be 
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started by selecting an initial value of a and updating 
it at every sample through the relationship, 
a (t)= a(t-l)+ A a (16) 

The control function is limited by the upper and lower 
limits and the pole shift factor should be such that it 
should be bounded by the reciprocal of the largest 
value of characteristic root of A (z _1 ). The latter 
requirement is satisfied by constraining the magnitude 
of a to unity. 

5. STATCOM Controller Performance 
for Single Machine System 

For the power system considered in figure 1, the input 
and output of the plant were considered to be the 
control of the modulation index (m) in the 
STATCOM circuit and the generator speed variation 
(Ago), respectively. In order to excite the plant, a 
sequence voltage steps and torque pulses in the 
regulator-exciter and generator shaft, respectively 
were applied. The diagonal elements of the initial 
covariance matrix P is assumed be 2xl0 5 ; the initial 
pole shift factor 0 and the forgetting factor of 1 were 
used. The starting values of all the parameters were 
considered to be 0.001 in all the simulations for 
consistency. The model order to be estimated was 
assumed to be 3. For solution of differential equations 
ode54 function from MATLAB 6.5 was used. The 
adaptive algorithm was developed by the authors. 




Figure 3. Generator speed deviation when excited with alternate 
torque steps 

Figure 3 shows the generator speed deviation with no 
control when excited by a sequence of toque steps of 
+5%, -5%, +5% and -5%. The nominal loading is 0.9 
pu at 0.9pf lagging. Figure 4 shows the variation of 
the generator speed with the pole-shift control applied 
to the identified process. 

From figures 3 and 4, it is apparent that the 
electromechanical transients are damped very well by 
the adaptive controller. The plant parameters are 
unknown at the start of the estimation process giving 
rise to very large overshoots. However, as the 
identification process progresses, the plant parameters 
are estimated more and more accurately to yield 
better updates of the pole shift factor, and hence 
providing better damping profiles. 




Figure 4 Speed deviation of the generator with adaptive pole- 
shift control 




Figure 5 The variation of the pole-shift factor with the 
progression of the adaptive process 




Figure 6 The variation of the a-parameters of the 
adaptively identified plant function 

Figure 5 shows the convergence of the pole shift 
factor as the estimation process progresses. The 
convergence of the {a} and {b} parameters of the 
estimated plant function are shown in figures 6 and 7, 
respectively. It can be observed that the estimation 
algorithm converges to the desired values rapidly. 
The convergence of the algorithm is independent of 
the initial choice of the pole shift factor a. 
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Figure 7 The variation of the b-parameters of the 
adaptively identified plant function 




Figure 9 Generator rotor angle variations following a three-phase 
fault for 0.1s with the proposed pole-shift control. The loadings are: 
a) 1 .2 pu, b) 0.9 pu, and c) 0.5 pu. 



6. Testing the STATCOM Controller 

A number of case studies were performed with the 
adapted system model and the pole shift parameters 
arrived at in the previous section. For a 50% input 
torque pulse on the generator for 0.1s, the rotor angle 
variations recorded for 3 operating conditions are 
shown in figure 9. These are for generator outputs of 
a) 1.2 pu , b) 0.9 pu , and c) 0.5pu. It can be observed 
that the damping properties are very good for the 
whole range of operation considered. Figure 10 shows 
the transient angle variations of the generator with the 
proposed adaptive control strategy for severe three- 
phase fault of 0.1s duration for the three loadings. It is 
to be noted that without control the system is under 
damped, in general, and unstable in some cases. 

7. Evaluation through PI Control 

The damping properties of the proposed coordinated 
robust controller were compared with a conventional 
PI controller in the STATCOM voltage magnitude 
control loop. The PI (or PID) controllers are normally 
installed in the feedback path. An additional washout 
is included in cascade with the controller to eliminate 
any unwanted signal in the steady state. The 
controller function in the feedback loop is written as, 




Figure 8 Generator rotor angle variations with the proposed 
pole-shift control for the loadings of a) 1.2 pu, b) 0.9 pu, and c) 0.5 
pu. The disturbance considered is a 50% torque pulse for 0.1s 



A pole-placement technique was used to determine 
the optimum gain settings (K P and K T ) of the 
controller. For a desired location of the dominant 
closed-loop eigen value X, the following equation is 
solved for K P and K I? 

H(X)= [Cftl-ApB]' 1 (18) 

H(X) is obtained from (17) for the desired X . 




Figure 10 PI controller configuration 



H(s) = [K P + 



Kj_ i iT„ 

s 1 + sT w 



(17) 




Figure 1 1 Generator rotor angle characteristics following a 
three-phase fault for 0.1s 
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Figure 11 shows the rotor angle variations of the 
synchronous generator following a three-phase fault 
with and without PI control for the 3 loading 
conditions considered in figures 8 and 9. The nominal 
loading is 1.01 pu. The values of K P and IQ, 
respectively are 5.957 and 10.892 for closed-loop 
system damping ratio of 0.3. The responses with 
dotted lines are with PI control. The figure shows 
comparison of the responses of the PI control with the 
pole-shift adaptive control of STATCOM voltage 
magnitude. While the PI control provides reasonable 
response for the nominal operating condition (curve 
b), for loadings farther away from nominal it is not 
satisfactory. The pole-shift adaptive control, on the 
other hand, provides very well damped stable 
response for all the 3 operating conditions considered. 

8. STATCOM Performance on a Multi- 
machine System 

A 4-machine power system with STATCOMs located 
at the middle of the transmission lines connecting 
each generator to the rest of the grid is shown in 
figure 12. The system dynamics involves a set of 
equations similar to (3) for each generator in addition 
to the network voltage-current relationships. The 
network quantities are transformed to individual 
machine d-q frames and a closed form state model for 
the multi-machine system is arrived at. 

The variable pole-shift adaptive control strategy was 
tested on the 4-machine power system given in figure 
12. The output was considered to be the change in 
power of generator #2 which is connected to the 
system through bus 9. The input signal is the 
modulation index of the STATCOM converter 
voltage. In order to excite the plant, a sequence of 
torque pulses was applied on the shaft of generator 
#2. 

Figure 13 shows the relative speed deviations of the 
generators with no control when excited by a 
sequence of toque steps of +5%, -5%, +5% and -5%. 
Figure 15 shows the variation of the generator speed 
with the adaptive variable pole-shift control applied to 
the identified process. The initial parameters of the 
adaptive algorithm were selected as in the single 
machine case. 

From figure 14, it is apparent that the 
electromechanical transients are damped very well by 
the adaptive controller. The plant parameters are 
unknown at the start of the estimation process giving 
rise to large overshoots. However, as the 
identification process progresses, the plant parameters 
are estimated more and more accurately to yield 
better updates of the pole shift factor, and hence 
providing better damping profiles. 

Figure 15 shows the convergence of the pole shift 
factor as the estimation process progresses. The 
convergence of the {a} parameters of the estimated 
plant function is shown in figures 16. It can be 
observed that the estimation algorithm converges to 
the desired values rapidly with minimum subsequent 



control expenditure. The convergence of the 
algorithm is independent of the initial choice of the 
pole shift factor a 




Figure 12 The multi-machine power system configuration 
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Figure 13 Generator speed deviation when excited with 
alternate torque steps 




Figure 14 Speed deviation of the generator with 
adaptive pole-shift control 
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Figure 15 The variation of the pole-shift factor with the 
progression of the adaptive process 

9. Testing the Multi-machine System 
Adaptive Controller 

A number of case studies were performed with the 
adapted system model and the pole shift parameters 
arrived at in the previous section. For a 50% input 
torque pulse on generator #2 for 0.1s, the rotor angle 
variations recorded for 2 operating conditions are 
shown in figure 17. The loadings are quite apart from 
each other. It can be observed that the damping 
properties are very good for the whole range of 
operation considered. It is to be noted that without 
control the system is very much under damped, and 
for this disturbance is on the verge of instability. 
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Figure 16 The variation of the a-parameters of the 
adaptively identified plant function 
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Figure 17 Relative rotor angle variations of the 
generators with the proposed pole-shift control for two 
loading conditions. The disturbance considered is a 100% 
torque pulse for 0.1s on generator 2. 

10. Conclusions 

An adaptive stabilizing control technique has been 
used to enhance the dynamic performance of a single 
machine as well as a multi-machine power system 
installed with synchronous static compensators. The 
control employed is the modulation index of one of 
the converter voltages. The proposed technique 
generates a stabilizing control on the basis of shifting 
the poles of the closed loop system towards the center 
of the unit circle in z-domain, thus providing more 
damping to the not so stable modes. The algorithm 
has been shown to converge to estimated parameter 
model rapidly for both the single machine as well as 
the multi-machine system. The on-line controller has 
demonstrated to provide very good damping to the 
electromechanical modes. The robustness of the 
control strategy was tested by simulating different 
types of disturbances including three phase faults 
covering a number of operating states. The adaptive 
strategy was also compared with optimized PI control 
response for a single machine system and was 
observed to provide much superior performance 
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